The Astrophysical Journal, 522:000-000, 1999 September 10 



COSMIC HISTORIES OF STARS, GAS, HEAVY ELEMENTS, AND DUST IN GALAXIES 

Yichuan C. Pei, S. Michael Fall, and Michael G. Hauser 

Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, Maryland 21218 

Received 1998 November 5; accepted 1999 April 21 
ABSTRACT 

We investigate a set of coupled equations that relates the stellar, gaseous, chemical, and radiation 
contents of the universe averaged over the whole population of galaxies. Using as input the available 
data from quasar absorption-line surveys, optical imaging and redshift surveys, and the COBE 
DIRBE and FIRAS extragalactic infrared background measurements, we obtain solutions for the 
cosmic histories of stars, interstellar gas, heavy elements, dust, and radiation from stars and dust 
in galaxies. Our solutions reproduce remarkably well a wide variety of observations that were 
not used as input. These include the integrated background light from galaxy counts from near- 
ultraviolet to near-infrared wavelengths, the rest-frame optical and near-infrared emissivities at 
various redshifts from surveys of galaxies, the mid-infrared and far-infrared emissivities of the local 
universe from the IRAS survey, the mean abundance of heavy elements at various epochs from 
surveys of damped Lyman-alpha systems, and the global star formation rates at several redshifts 
from Ha, mid-infrared, and submillimcter observations. The chemical enrichment history of the 
intergalactic medium implied by our models is also consistent with the observed mean metal content 
of the Lyman-alpha forest at high redshifts. We infer that the dust associated with star forming 
regions is highly inhomogeneous and absorbs a significant fraction of the starlight, with only 41—46% 
of the total in the extragalactic optical background and the remaining 59 — 54% reprocessed by dust 
into the infrared background. The solutions presented here provide an intriguing picture of the 
cosmic mean history of galaxies over much of the Hubble time. In particular, the process of galaxy 
formation appears to have undergone an early period of substantial inflow to assemble interstellar 
gas at z J> 3, a subsequent period of intense star formation and chemical enrichment at 1 <J z < 3, 
and a recent period of decline in the gas content, star formation rate, optical stellar emissivity, and 
infrared dust emissivity at z <J 1. 

Subject headings: cosmology: diffuse radiation — galaxies: evolution — quasars: absorption lines 
— infrared: galaxies 



1. INTRODUCTION 

In the past few years, galaxies have been observed to 
redshifts approaching five, when the universe was less 
than 10% of its present age. A goal of several recent 
studies has been to determine the average evolution of 
the stellar and interstellar contents of the entire popula- 
tion of galaxies. Individual galaxies, spanning a wide 
range of sizes, masses, and morphological types, cer- 
tainly had different, possibly complicated histories of 
star formation, gas consumption, and chemical enrich- 
ment. Studying the average properties of galaxies has 



the virtue of simplicity and is relevant to several im- 
portant problems, such as the extragalactic background 
radiation and the chemical enrichment of the intergalac- 
tic medium. Quasar absorption-line and optical imaging 
and redshift surveys of galaxies, when combined with 
cosmic infrared background measurements, promise to 
reveal the unbiased global history of star formation in 
galaxies. This is the main objective of this paper. 

Quasar absorption lines have been used for more than 
a decade to probe the interstellar content of galaxies at 
high redshifts. The damped Lya systems, a subset of the 
absorption-line systems with highest HI column densi- 
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ties, are of particular interest because they dominate the 
mass density of neutral gas in the universe and probably 
represent the progenitors of present-day galaxies (Wolfe 
et al. 1986), although the exact nature of these objects 
(i.e., their morphologies, sizes, and space density) is a 
topic of current speculation and investigation. Surveys 
of damped Lya systems along many lines of sight provide 
information on the mean comoving density of neutral 
hydrogen and the mean abundances of heavy elements, 
dust, and molecules at various redshifts. The picture 
that emerges, while still incomplete, already reveals sev- 
eral important results. First, the damped Lya systems 
at z ~ 3 contain enough cool, neutral gas to form all 
of the stars visible today (Wolfe et al. 1995). Second, 
the observed decrease of the gas content from z « 3 
to z = plausibly reflects the conversion of gas into 
stars (Lanzctta, Wolfe, & Turnshek 1995). Third, the 
damped Lya systems at z « 3 are still in their infancy 
of chemical enrichment, with low metallicity (Pettini et 
al. 1994), low dust-to-gas ratio (Pei, Fall, & Bechtold 
1991) and few molecules (Levshakov et al. 1992). These 
observations, when tied together by models of cosmic 
chemical evolution, predict high rates of star formation 
at 1 <; z <^ 2 and a rapid decline toward z = (Pei & 
Fall 1995). 

Recently, optical imaging and redshift surveys have 
given us a direct view of the stellar content of the dis- 
tant universe (Lilly et al. 1995; Steidel et al. 1996; Ellis 
et al. 1996; Cowie et al. 1996; Williams et al. 1996). 
An important quantity that can be measured from these 
surveys is the mean rest-frame ultraviolet emissivity as a 
function of redshift, which potentially traces the global 
history of star formation in the universe. The optical 
data suggest that the average rate of star formation in- 
creases with time at z ^ 2 (Madau et al. 1996), peaks 
near z « 1 — 2 (Connolly et al. 1997), and then declines 
sharply to z = (Lilly et al. 1996), although the ap- 
parent increase from z ~ 4 to z m 3 is quite uncertain 
(Steidel et al. 1999). This picture is broadly consistent 
with the history of star formation inferred from models 
of cosmic chemical evolution with input from surveys of 
damped Lya systems (Fall 1998). Thus, cosmic chem- 
ical evolution provides a key link between the stellar, 
gaseous, and chemical constituents of galaxies. 

All-sky surveys from the DIRBE and FIRAS instru- 
ments on the COBE satellite have provided us with an- 
other view of the integrated emission history of the uni- 
verse from far-infrared to millimeter wavelengths (Puget 
et al. 1996; Schlegel, Finkbeiner, & Davis 1998; Hauser 
et al. 1998; Fixsen et al. 1998). The recent DIRBE and 
FIRAS detections of the cosmic infrared background, 
when compared with the Hubble Deep Field survey, in- 
dicate that a significant fraction of the radiation released 



by stars has not been seen in the ultraviolet and opti- 
cal region of the extragalactic background, but rather at 
far-infrared wavelengths (Hauser et al. 1998). The large 
amount of far-infrared radiation seen by the DIRBE and 
FIRAS appears to imply that the global star formation 
rate currently deduced from optical surveys is severely 
biased (Dwek et al. 1998; Blain et al. 1999; Calzetti & 
Hcckman 1999). This bias is not entirely surprising be- 
cause galaxies undergoing bursts of star formation are 
known not only to be sources of intense ultraviolet radi- 
ation but also of far-infrared radiation. The cosmic in- 
frared background records the cumulative infrared emis- 
sion from galaxies at all redshifts and therefore provides 
an important constraint on the global history of star 
formation. 

In this paper, we combine the available data from 
quasar absorption-line surveys, optical imaging and red- 
shift surveys, and extragalactic background measure- 
ments to derive the cosmic histories of star formation, 
gas consumption, and chemical enrichment in galaxies. 
Our approach is an extension of the studies of Pei & Fall 
(1995) and Fall, Chariot, & Pei (1996). These are based 
on a set of equations that link the comoving densities 
of stars, interstellar gas, heavy elements, and dust, av- 
eraged over the whole population of galaxies. The data 
from quasar absorption-line surveys and optical imaging 
and redshift surveys are taken as inputs, and the effects 
of dust are self-consistently coupled to the production 
of heavy elements. A key constraint in the present pa- 
per is imposed by the recent DIRBE and FIRAS mea- 
surements of the cosmic infrared background. With a 
self-consistent treatment of the absorption and reradia- 
tion of starlight by dust, the problem is then to solve 
the set of coupled chemical equations so that the rera- 
diated starlight by dust matches the background mea- 
surements. Section 2 presents the basic method, states 
the key assumptions, and identifies the necessary obser- 
vational inputs. In § 3, we derive the global histories 
of stars, interstellar gas, heavy elements, dust, radiation 
from stars and dust, and baryon flow between the in- 
terstellar and intergalactic media. In § 4, we discuss the 
major uncertainties in the models and some of the impli- 
cations of our results. Our conclusions are summarized 
in § 5. Throughout the paper, we adopt a cosmology 
with A = 0, q = 1/2, and H = 50 km s _1 Mpc -1 . 

2. MODELS 

We are interested in the global evolution of stars, in- 
terstellar gas, heavy elements, dust, and radiation from 
stars and dust, all averaged in comoving volumes large 
enough to contain many galaxies. The basic formalism 
that can be used to describe this evolution has been de- 
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veloped by Lanzetta et al. (1995), Pci & Fall (1995), and 
Fall et al. (1996). Our previous work showed that the 
emission history of galaxies could be inferred from their 
absorption history with the aid of stellar population syn- 
thesis models, and that self-consistent solutions of the 
equations of cosmic chemical evolution predicted high 
rates of star formation at 1 <; z <^ 2 and declining rates 
at ^ z ^ 1 • Here, we adopt nearly the same formalism 
in describing cosmic chemical evolution but without any 
assumption about the baryon flow rate between the in- 
terstellar and intergalactic media, and develop some new 
formulae for the absorption and reradiation of starlight 
by dust. With observations now available on the ab- 
sorption and emission histories of galaxies, we can use 
both of these as inputs, while the reradiation of starlight 
absorbed by dust is constrained by the measured cosmic 
infrared background radiation. In this way, we obtain 
the unbiased absorption and emission histories of galax- 
ies as well as the evolution of baryons within galaxies. 

2.1. Cosmic Chemical Evolution 

The equations of cosmic chemical evolution govern 
the global rates of star formation, gas consumption, and 
chemical enrichment. For a population of galaxies, these 
rates can be expressed most simply in terms of the mean 
comoving densities of stars, interstellar gas, and ele- 
ments heavier than He in the interstellar medium, Q Sl 
il g , and f2 m , all measured in units of the present criti- 
cal density (p c = 3Hq/8itG), and the mean abundance 
of heavy elements in the interstellar media of galaxies 
Z = Q m /fl g . As defined here, f2 s , fl gi f2 m , and Z all 
pertain to galaxies themselves, not to the intergalactic 
medium. On the timescale of interest here (the Hubble 
time), any delayed recycling of stellar material to the 
interstellar medium can be neglected. We can therefore 
write, in the approximation of instantaneous recycling 
(and Z«l), 

n g + n s = h f , (i) 

Sl g Z-y(l a = (Z f -Z)(l f , (2) 

where y is an "effective yield" of heavy elements and the 
dots denote differentiation with respect to proper time. 
The source terms on the right-hand sides of equations 
(1) and (2), with subscript /, represent the inflow or 
outflow of gas with mean metallicity Zf at a net comov- 
ing rate (if to or from galaxies. These equations can 
be derived by summing the corresponding equations for 
the chemical evolution of all the individual galaxies in a 
large comoving volume (see Tinsley 1980). 

The definition of the effective yield requires some care. 
It could be defined as the mean nucleosynthetic yield 
of different galaxies weighted by their star formation 



rates. In this case, however, equation (2) would only be 
valid for a population of chemically homogeneous galax- 
ies (i.e., all galaxies at a given redshift have the same 
metallicity) . The reason for this is that the average over 
galaxies leading to equation (2) contains a residual term 
due to the metallicity spread of different galaxies, which 
can simply be absorbed into the definition of y. In this 
case, equation (2) is also valid for a population of chem- 
ically inhomogeneous galaxies. This definition of the 
effective yield is entirely appropriate in our models, be- 
cause we regard y as an output parameter fixed by the 
mean metallicity at the present epoch (see below). 

We now consider a simple case of the inflow or out- 
flow of gas between galaxies and their environment. In- 
flow (Clf > 0) allows for the addition of unprocessed 
material (Zf = 0) into galaxies from the intergalactic 
medium, while outflow (fi f < 0) allows for the removal 
of processed material (Zf = Z) from galaxies by super- 
nova explosions, galactic winds, collisions, and so forth. 
Thus, we write 



z f = z[i-e((if)}, 



(3) 



where theta is the unit step function, defined as 9(x) = 1 
for x > and 9(x) = for x < 0. If a significant 
outflow occurs, it can enrich the intergalactic medium 
with heavy elments. However, since the gas reser- 
voir in the intergalactic medium is sufficiently large, its 
overall metallicity is likely to remain lower than that 
of the interstellar medium within galaxies. Therefore, 
the assumption of only metal-free inflow should remain 
roughly valid. In principle, there is a possibility that 
both inflow and outflow could occur at the same time 
and with the same rate so that the net flow is zero, in 
which case equation (3) would be invalid. In practice, 
as long as the net flow is dominated by cither inflow or 
outflow, equation (3) should be a good approximation. 

Equations (1) and (2) can be solved by specifying the 
initial conditions at t = or the boundary conditions at 
the present epoch, z = 0. We assume that, at some early 
time, before galaxies form, there are no stars, interstellar 
gas, or heavy elements, i.e., fl s = il g = Cl m = at 
t = 0. With these initial conditions and equation (3), 
the solutions to equations (1) and (2) are 
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The formation and early evolution of galaxies naturally 
requires an initial inflow, (if > 0. Any time variation in 
y is likely to be weaker than that in (l s , which varies by 
more than an order of magnitude over <; z <; 5. As 
an approximation, we assume that the effective yield y 
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is constant in time. Since Z and y always appear in the 
combination Z/y, we can fix y by imposing a boundary 
condition on the mean metallicity Z at z = 0. The quan- 
tity Z(0), which refers to an average over all galaxies at 
the present epoch, is not known precisely. We simply 
assume that the metallicity in the solar neighborhood is 
typical and adopt Z(0) = Z©. The values of Z(0) and 
all other input parameters in our models are listed in 
Table 1. 

Our choice of the initial conditions for Sl s and Sl g has 
a special meaning in the models of cosmic chemical evo- 
lution: the quantity Slf at any redshift is always the sum 
of the comoving densities of stars and interstellar gas, 
and hence represents the comoving density of baryons 
within galaxies. The evolution of the baryon content 
in galaxies is regulated by the star formation rate and 
the gas consumption and inflow or outflow rates. If the 
star formation rates were relatively small at some early 
epochs, i.e., Cl s <C Cl g , the evolution of galaxies would 
be dominated by the build up of interstellar gas, i.e., 
Clf s=s Clg . On the other hand, if the star formation rates 
were relatively high, this would dominate the evolution. 
In fact, the process of galaxy formation might proceed 
through two phases: (1) the accretion of baryons during 
an early period of gas assembly, and (2) the production 
of heavy elements and consumption of gas during a sub- 
sequent period of star formation. 

The total abundance of heavy elements Z in equation 
(5) includes metals in both the gas and solid phases of 
the interstellar medium. The production of dust is there- 
fore included naturally in our approach. We assume for 
simplicity that a fixed fraction of the heavy elements in 
the interstellar medium is locked up in dust grains. In 
this case, the comoving density of dust fid is related to 
the comoving density of heavy elements Sl m by 



fid — dm,flr< 



dm ZSl 
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(6) 



with a constant dust-to- metals mass ratio d m . In ad- 
dition to some theoretical expectations for a constant 
d m (Dwck 1998), there is also empirical evidence for 
this in nearby galaxies and the damped Lya systems 
at high redshifts. For the Milky Way and the Large and 
Small Magellanic Clouds, the dust-to-metals ratio is 0.51 
(GAL), 0.36 (LMC), and 0.46 (SMC) (Pei 1992; Luck & 
Lambert 1992), roughly constant even though the dust- 
to-gas ratios and metallicities in these galaxies vary by 
nearly an order of magnitude. Using the measurements 
of the gas-phase depletion of Cr relative to Zn in a sam- 
ple of damped Lya systems by Pettini ct al. (1997a), we 
derive mean values d m = 0.45 at 0.7 < z < 1.0, 0.37 
at 1.0 < z < 2.0, and 0.44 at 2.0 < z < 2.8, with the 
assumptions that the intrinsic Cr/Zn is solar, that Cr 
is depleted onto grains, while Zn is not, and that the 



dust has Galatic chemical composition and SMC-type 
extinction curve. These results are consistent with es- 
timates based on the observed gas-phase depletions of 
other elements in damped Lya systems (Kulkarni, Fall, 
& Truran 1997; Vladilo 1998). Figure 1 shows our es- 
timates of d m as a function of redshift, where the data 
point at z — is the mean for the Milky Way, LMC, 
and SMC. Evidently, the mean dust-to-metals ratio is 
roughly constant at d m — 0.45 over a significant range 
of redshifts, < z < 3. 




FIG. 1. — Mean dust-to-metals mass ratio d m as a 
function of redshift. The data point at z = is the mean 
in the Milky Way, LMC, and SMC (Pei 1992; Luck & 
Lambert 1992), while the others are derived from mea- 
surements of the gas-phase depletion of Cr relative to 
Zn in damped Lya systems (Pettini et al. 1997a). The 
error bars represent la uncertainties in the means. 



The equations presented above for cosmic chemical 
evolution are quite general and can be used to study the 
main baryonic ingredients of galaxies and their evolu- 
tion. The global quantities of interest here include the 
comoving densities of stars (Sl s ), interstellar gas (Sl g ), 
heavy elements (Sl m ), dust (Sid), an d baryons (Slf) in 
galaxies, and the comoving rates of star formation (Cl s ), 
gas accretion (Cl g ), metal enrichment (Cl m ), dust pro- 
duction (Sid), and baryon flow (Slf), all as functions 
of redshift. With equations (4)- (6) and straightforward 
time differentiation and integration between the comov- 
ing density and rate pairs, there are only two (func- 
tional) degrees of freedom in the problem, namely, in- 
puts of any two independent quantities as functions of 



4 



rcdshift allow all others to be derived. We choose £l g (z) 
and Sl s (z) as the two inputs, since both can be traced 
over a wide range of rcdshifts by observations. 

2.2. Absorption History 

The damped Lya systems have HI column densities in 
excess of 2 x 10 20 cm -2 and account for more than 80% 
of the neutral gas in the universe (Lanzetta et al. 1991). 
The other quasar absorption-line systems, those with 
HI column densities less than 10 20 cm~ 2 , probably con- 
tain more gas in total than the damped Lya systems, 
but this gas is diffuse and highly ionized. The damped 
Lya systems are usually interpreted as the ancestors of 
present-day galaxies, but it is not yet clear whether they 
are mainly the progenitors of galactic disks or spheroids, 
or dwarf galaxies. However, the comoving density of HI, 
the quantity of interest here, can be computed directly 
from the statistics of the absorption along random lines 
of sight and is therefore independent of the morpholo- 
gies, sizes, and other properties of the damped Lya sys- 
tems. The main restriction on our results is that they do 
not include any galaxies that consumed their neutral gas 
before z w 4, the highest redshift probed systematically 
by existing surveys of damped Lya systems. 

2.2.1. Comoving Density of Gas 

The interstellar content of galaxies consists mainly of 
hydrogen and helium (in neutral atomic, ionized, and 
molecular forms); thus 

n g = /uQ H = M^hi + ftffll + ^h 2 ), (7) 

where /i is the mean particle mass per H atom, and HI, 
HII, and H2 denote neutral atomic, ionized, and molec- 
ular hydrogen, respectively. The damped Lya systems 
are probably mainly neutral, because their HI column 
density exceeds N ^ 2 x 10 20 cm~ 2 and their metal ab- 
sorption lines usually show low ionization states. They 
also appear to have low abundances of molecules, at least 
at high rcdshifts (Levshakov et al. 1992; Ge & Bechtold 
1997). As an approximation, we ignore HII and H 2 but 
include He: 

Q g = /xO H i, (8) 

with /1 = 1.3, appropriate for 23% primordial He by 
mass from Big Bang nucleosynthesis (Hogan, Oliver, & 
Scully 1997). While the neglect of H 2 is probably rea- 
sonable at high rcdshifts, it is less accurate toward z — 
0. In the local universe, the observed H 2 and HI masses 
of 154 late-type galaxies (Sa to Im) give an averaged 
H2-to-HI mass ratio of roughly 0.8, implying 45% of H 
in H 2 (Young & Scoville 1991). Thus, it appears that 
equation (8) is approximately valid (to within a factor 



of two) over much or all of cosmic history. 

The comoving density of HI as a function of redshift 
is given by 

nm(s) = dNN f(N, z), (9) 

where ijih is the mass of the H atom, and / is the 
distribution of HI column densities, defined such that 
f(N, z)dNdX is the number of absorbers along a ran- 
dom line of sight with HI column densities between N 
and N + dN and redshift paths between X and X + dX 
[with dX = (l + z)(l + 2q a z)- 1 / 2 dz for A = 0] (Lanzetta 
et al. 1991). The existing samples of damped Lya sys- 
tems are all derived from optically selected quasars. As 
a result of the obscuration caused by dust within the 
damped Lya systems, these samples are biased to some 
degree (Fall & Pei 1993). Thus, we introduce the ob- 
served comoving density of HI 

*W*) = r dNN f {N, z), (10) 

where f a represents the observed distribution of HI col- 
umn densities in damped Lya systems. In the limit of 
an infinitely large sample of absorbers, the true and ob- 
served distributions, / and f a , sample all impact param- 
eters and orientations of galaxies, and the resulting co- 
moving densities, Ohi and f^Hio, can be derived without 
any knowledge of the sizes or shapes of galaxies. 

Both f^Hi and flnio are averages over the population 
of damped Lya absorbers, but the former pertains to 
random lines of sight, while the latter pertains to the 
lines of sight to optically selected quasars. Since the 
missing absorbers in an optically selected sample tend 
to be those with the highest N, there can be large dif- 
ferences between fini and Hhio and this can affect the 
resulting picture of cosmic chemical evolution (Pei & 
Fall 1995). We neglect another potential bias, that due 
to the magnification of quasars by gravitational lenses 
associated with foreground damped Lya systems, be- 
cause this appears to be a relatively minor effect in the 
samples from which finio has been estimated (Le Brun et 
al. 1997; Pernact al. 1997; Smcttc ct al. 1997). Since the 
dust-to-gas ratios are not known individually for many 
damped Lya systems, the conversion from fimo to fim 
must be made statistically. For this purpose, we intro- 
duce a mean correction factor Q, defined by the relation 

n m ( z ) = n mo (z)Q(z). (11) 

2.2.2. Missing Absorbers 

The relation between the true and observed distribu- 
tions of HI column densities can be derived under two 
idealizations: that the absorbers are obtained from a 
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magnitude limited sample of bright quasars, and that 
all absorbers at a given redshift have the same dust-to- 
gas ratio. In this case, we have (Fall & Pci 1993) 

f(N,z) = f {N,z)e X p\pT(N,z)], (12) 

t(N, z) = k m (z)m H N Ke [\ e /(l + z)}, (13) 

k m (z) = n d (z)/n m (z), (14) 

where (3 is the power-law index of the bright part of 
the quasar luminosity function [i.e., 4>(L,z) oc L^ 13 ^ 1 ], 
t is the extinction optical depth at z = 0, k m is the 
mean dust-to-HI mass ratio, K e is the mass extinction 
coefficient evaluated at a wavelength A = A e /(1 + z), 
and A e is the effective wavelength of the limiting mag- 
nitude. Equations (12)-(14) are also valid for combi- 
nations of samples with different limiting magnitudes, 
provided these are brighter than V <^ 20, and they are 
good approximations when there is a dispersion in the 
dust-to-HI ratio, provided this satisfies cr(log£; m ) <J 1. 
Both conditions appear to be true for the existing sam- 
ples of damped Lya systems and background quasars 
(see § 2.2.3 and § 2.2.4). 

To derive the correction factor Q, we must adopt a 
specific form for the observed distribution of HI column 
densities. A common approximation is a power law, 
fo{N,z) oc iV-T, with 7 « 1.7 (Lanzetta et al. 1991). 
Unfortunately, this gives a divergent total mass in equa- 
tion (10), unless an upper cutoff in TV is imposed. Thus, 
following Pei & Fall (1995), we adopt the gamma distri- 
bution 

f a (N, z) = (f*/N*)(N/N*)^ expi-N/N.), (15) 

where /* and N* specify the normalization and "knee" of 
the function. This parametric form provides a solution 
to the problem of diverging mass and a better fit to 
the data than the power-law form (Storrie-Lombardi et 
al. 1996a). In this case, the true distribution / from 
equation (12) is also a gamma distribution 

f(N,z) = (f t /N t )(N/N t )-~<eM-N/N t ), (16) 

with f t = Q^.U and N t = Q^N*. For 7 < 2, 
the observed and true comoving densities of HI are sim- 
ply fi fflo = (8iGm H /3cff )r(2 - 7 )/*iV, and Sl m = 
(87rGm H /3c_ff )r(2 — j)f t N t . Combining equations (9)- 
(15) and eliminating in favor of 7 and /*, we obtain 

Q = (ir, +Q^, (17) 

T * (Ae ' Z) = 8^Gf~^T) Qd{z)Ke[Xe/{1 + Z)l (18) 
where V denotes the standard gamma function. Equa- 
tion (17) is an implicit relation from which Q can be 
obtained as a function of /3r* and 7. For 7 = 1, we 
have Q = 1 + /3r* and, for 7 = 3/2, we have Q = 



[l + (/3r*/2) 2 ] 1/2 + (/3r*/2). 

2.2.3. Inputs from Damped Lya Surveys 

The formulae presented above enable us to compute 
the true Q, g from the observed comoving density of HI, 
finio- This relation is 

Q g (z) = fjQmo(z)Q(z), (19) 

where (i accounts for He and Q accounts for the missing 
damped Lya systems in an optically selected sample of 
quasars. The correction factor Q depends on the mean 
comoving density of dust Qd, the extinction coefficient 
of grains K e , the bright-end slope of the quasar lumi- 
nosity function f3, the effective wavelength of the quasar 
sample A e , the slope of the absorber column density dis- 
tribution 7, and the characteristic sky coverage of the 
absorbers /*. Because Q is linked to fid, our solutions 
of the equations of cosmic chemical evolution are ob- 
tained itcrativcly (see § 3). In the following, we adopt 
^Hio(z) at various rcdshifts from surveys of the damped 
Lya systems and estimate the values of the parameters 
(/?, A e , 7, /*) appropriate for these surveys. 

The observed HI content in the damped Lya systems 
has been traced over the range of redshift 0.008 < z < 
4.7, in surveys by Wolfe et al. (1986, 1995), Lanzetta 
et al. (1991, 1995), and Storrie-Lombardi et al. (1996b). 
We adopt the most recent determinations of Qyu (z), 
derived by Storrie-Lombardi et al. (1996b) assuming a 
particular cosmology with A = 0, q = 1/2, and Ho = 
50 km s _1 Mpc -1 . These estimates, listed in Table 2 
with their errors, are based on 44 damped Lya systems 
obtained from a combined sample of 366 optically se- 
lected quasars. The statistical uncertainty in S]hi c (z) is 
roughly 50%. The present value Ohi o (0) is not directly 
determined from absorption in damped Lya systems, 
and, as a result, is not available as input to our models. 
Instead, we appeal to observations of 21 cm emission for 
the HI content of galaxies in the local universe (§ 2.2.6). 

The existing samples of damped Lya systems have 
been obtained from samples of quasars selected by V < 
18.5 (Wolfe et al. 1986; Lanzetta et al. 1991), B < 18.75 
(Wolfe et al. 1995), V £ 17 (Lanzetta et al. 1995), 
and R < 19.5 (Storrie-Lombardi et al. 1996a). Thus, 
they do not conform to the idealization of a single lim- 
iting magnitude under which equation (12) was derived. 
Fortunately, all of the samples are brighter than V w 
20. In this case, the luminosity function of quasars is 
well approximated by a single power law (Hartwick & 
Schadc 1990; Pei 1995). The limiting magnitudes can- 
cel out and do not appear explicitly in equation (12), 
except through the effective wavelength A e . Therefore, 
our formulae for Q are also valid for combinations of 
samples with different limiting magnitudes. The value 
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of (3 can be determined most simply from the differential 
counts A(m) of quasars versus apparent magnitude m, 
because 4>{L,z) oc _L~' 3 ~ 1 implies A{m) oc l0°- 4 ' 3m . Fig- 
ure 2 shows the observed _B-band counts from Hartwick 

6 Schade (1990). The line is a weighted, least-squares 
fit to the data at B < 20, with the x 2_n tt e( i value of 
(3 listed in Table 1. At z < 3.5, most of the absorbers 
come from the ^-selected quasars, while at z > 3.5, all 
of them come from the i?-selected quasars. Accordingly, 
we adopt A e = 5500 A at z < 3.5 and A e = 6500 A at 
z > 3.5. The value of A e at z = docs not enter our 
calculations but implies a specific value of f2m o (0) for 
the purpose of comparison (§ 2.2.6). 

The current data for the observed distribution f a in 
the damped Lya systems are sparse but nevertheless 
indicate that the evolution of f a is strongest for the 
high-A systems (Wolfe et al. 1995; Storrie-Lombardi et 
al. 1996a). In the gamma distribution model, this is 
equivalent to varying ./V* with respect to z, and hence the 
evolution of SIhio is predominantly due to A* rather than 

7 and /*. Thus, it is a good approximation to regard 
7 and /* as constants. We determine these parameters 
by fitting the cumulative form of / , i.e., the observed 
number n Q (> A) of absorbers with HI column densi- 
ties greater than A. The results are shown in Figure 3, 
where the data (in the stepped line and the circle) are 



from Storrie-Lombardi et al. (1996a). The curve is the 
gamma distribution (integrated over the redshift path 
of all absorbers in the sample), with the fitted values 
of 7 and /* listed in Table 1. The circle at A = 1.6 x 
10 17 cm~ 2 is the number of Lyman-limit systems that 
would be detected along the same redshift path as in the 
damped Lya sample. We have included it in our fits to 
reduce further the uncertainty in 7 and /* , especially the 
former. The values of 7 and /* derived here are similar to 
those obtained by Storrie-Lombardi et al. (1996a) from 
a maximum likelihood fit. 

2.2.4- Dispersion in the Dust-to-Gas Ratio 

The main idealization in our derivation of the mean 
correction Q is the assumption that all absorbers at a 
given redshift have the same dust-to-gas ratio. This 
gives a simple relation between the true and observed 
distributions in equation (12) so that Q can be evalu- 
ated using f (or /). It also links the dust-to- HI ratio 
k m to the comoving density of dust in equation (14) 
so that Q is coupled self-consistcntly to the equations 
of chemical evolution. The damped Lya systems do 
show a significant dispersion in their dust-to-gas ratios, 
as indicated by the observed Cr/Zn ratios (Pettini et 
al. 1997a). A dispersion in the dust-to-gas ratio would 
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FIG. 2. — Differential counts of quasars per square de- 
gree versus B magnitude. The data points were derived 
by Hartwick & Schade (1990) from a combined sample 
of more than 1000 quasars over < z < 3.3. The line is 
a weighted, least-squares fit to the observations at B < 
20, relevant to the samples of quasars from which the 
damped Lya systems were selected. 
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FIG. 3. — Cumulative number of absorption systems 
with HI column densities greater than A. The stepped 
line represents the data from existing surveys of damped 
Lya systems (with a total redshift path AX = 421), 
while the point is the mean number of Lyman limit sys- 
tems that would be detected along the same redshift 
path (Storrie-Lombardi et al. 1996b). The curve shows 
the fit of the cumulative form of the gamma distribution. 
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introduce additional terms of N in the relation between 
/ and f , rather than the simple relation in equation 
(12) (see the Appendix of Fall & Pei 1993). Nearby 
galaxies have cr(log/c m ) « 0.4 — 0.5 (Pei 1992) and the 
observed damped Lya systems have <r(logfe m ) « 0.6 at 
0.7 < z < 2 and 0.6 at 2 < z < 2.8 (Pcttini et al. 1997a). 
If the true dispersion were not much larger than the ob- 
served one, the additional terms in the relation between 
/ and f would be negligible (less than roughly 20% at 
N = 10 21 cm" 2 for k m = 0.1 Z Q ). Thus, the mean 
correction Q derived here is not sensitive to a modest 
dispersion in the dust-to-gas ratio [i.e., cr(logfc m ) <; 1]. 

A dispersion in the dust-to-gas ratio could, however, 
affect the observed mean dust-to-gas ratio and metal- 
licity in damped Lya absorbers. The observed mean 
dust-to-gas ratio could be underestimated, because the 
absorbers with highest dust-to-gas ratios at fixed HI col- 
umn density would be preferentially missing from opti- 
cally selected quasars. Using the formulae of Fall & Pei 
(1993), we estimate that the observed mean dust-to-gas 
ratio would be 30% — 50% lower than the true mean for 
cr(logfc m ) w 0.4 - 0.6 (at k m = 0.1 Z Q ). If the heavy 
elements in individual galaxies were correlated with the 
dust, the observed mean metallicity would also be un- 
derestimated by a similar amount. These biases, while 
not negligible, are comparable to the observational un- 
certainties in the derived mean dust-to-gas ratios and 
metallicities in the current samples of damped Lya sys- 
tems. We emphasize that our models of cosmic chemi- 
cal evolution do not use as input the observed comoving 
densities of dust and metals, and consequently they are 
not affected by a dispersion in the dust-to-gas ratio. The 
only complication is that, since no corrections are made, 
this could affect the comparison of our model outputs 
with the observed mean dust-to-gas ratio and metallic- 
ity in the damped Lya systems. 

2.2.5. Optical Properties of Grains 

The optical properties of interstellar dust grains at 
high redshifts are largely unknown. We must therefore 
appeal to the properties of galaxies at the present epoch. 
The Milky Way, LMC, and SMC are three nearby galax- 
ies with well observed interstellar extinction curves. To 
derive self-consistently the wavelength dependence of the 
albedo (relevant for the absorption of starlight by dust 
discussed in § 2.3), we adopt the Draine & Lee (1984) 
grain model but with the proportions of graphite and 
silicates adjusted so as to fit the mean empirical ex- 
tinction curves in the Milky Way, LMC, and SMC (Pei 
1992). Figure 4 shows the resulting opacities, n e (mass 
extinction coefficient in the upper panel) and n a (mass 
absorption coefficient in the lower panel), as functions 
of wavelength in the three galaxies. The range in n e be- 



tween the Milky Way, LMC, and SMC is roughly 30%, 
whereas the range in n a can be as large as a factor of 
two, especially at optical wavelengths. The LMC and 
SMC have lower abundances of metals and dust than the 
Milky Way and may be more representative of galaxies 
at high redshifts. There is some evidence for this from 
the weak or absent 2200 A "absorption" feature in the 
damped Lya systems (Pei et al. 1991). In our calcu- 
lations, we adopt the LMC-type dust, representing an 
intermediate case between the Galactic-type and SMC- 
type dust. We will later discuss the effects of different 
n a on our results (§ 4.2). 




A/ fim 

FIG. 4. — Extinction and absorption coefficients (per 
unit mass) of grains in the Milky Way, LMC, and SMC. 
These are based on the Draine & Lee (1984) model but 
with the proportions of graphite and silicates adjusted 
to fit the observed mean extinction curves in the three 
galaxies (Pei 1992). 

2.2.6. Neutral Gas in the Local Universe 

The existing samples of damped Lya systems con- 
tain only a few absorbers at z <; 2, causing nearly 60% 
uncertainty in the corresponding estimates of Q,m.o{z) 
in Table 2. The value of Slnio at z = 0, if extrapo- 
lated from these entries, is even more uncertain. To 
minimize this uncertainty, we note that averages over 
the HI content of galaxies in the local universe give 
Q HI (0) = (3.8 ± 0.9) x 10- 4 (Fall & Pei 1993; Rao 
k Briggs 1993; Zwaan et al. 1997). This estimate of 
Ohi(0) is based on 21 cm emission from galaxies, and 
unlike fiHio(O), it is not affected by dust obscuration. 
Thus, we adopt this to fix the present true comoving 
density of gas in our models, listed in Table 2 as the 
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first entry under Q, g (z). Galaxies at z = contain 
dust and hence cause some difference between f2Hi(0) 
and f2Hio(0). The value of f2m(0) adopted here implies 
^Hio(O) = (2.2 — 2.7) x 10~ 4 in our models, consistent 
with the range of n mo (0) = (2.1 - 3.0) x 1(T 4 that 
would be extrapolated from the data in the damped Lya 
systems [with a form log^Hio^) a z or logf^Hio^) cx 
log(l + z) at z £2]. 



2.3. Emission History 

We now turn to the second input to our models, 
Cl s (z). Deep optical imaging surveys probe the rest- 
frame ultraviolet radiation of high-redshift galaxies and 
the rest-frame optical radiation of low-redshift galax- 
ies. When combined with redshift estimates, they trace 
the stellar emission history of galaxies. The ultraviolet 
light is produced mainly by short-lived, massive (O and 
early B) stars and can be used to infer the instantaneous 
star formation rate, nearly independent of the prior star 
formation history. However, ultraviolet radiation may 
be strongly absorbed by dust, and the star formation 
rate may be severely underestimated if no corrections 
are made. Since the absorbed starlight is reradiated by 
dust, the recent measurements of the cosmic infrared 
background by the DIRBE and FIRAS provide a crucial 
constraint, which allows us to correct for any radiative 
energy that is missing in the optical data. The infrared 
background is the cumulative result of emission at all 
redshifts. When combined with the chemical equations, 
the redshift information is then decoded from the inte- 
grated background constraint. This section provides a 
self-consistent treatment of the absorption and reradia- 
tion by dust in order to obtain unbiased estimates of the 
global rates of star formation. 

2.3.1. Emissivity and Background Intensity 

We are interested in the global evolution of the radia- 
tion from stars and dust in galaxies. This evolution can 
be described naturally in terms of the comoving emis- 
sivity E v {z) as a function of redshift, defined here as the 
power radiated per unit frequency per unit comoving 
volume at the redshift z. When computing E Ul we in- 
clude only photons that have escaped from galaxies and 
neglect the absorption of infrared photons. This emissiv- 
ity is then a sum of the direct (but attenuated) starlight 
and the reradiated starlight from dust in galaxies: 



E v = (1 — A V )E SU + Ed v , 



(20) 



where A v is the mean fraction of photons absorbed by 
dust in the interstellar medium, and E sv and E^v are the 
stellar and dust emissivities, respectively. It is clear that 
E sv dominates at ultraviolet, optical, and near-infrared 



wavelengths, while E& v dominates at mid-infrared, far- 
infrared, and submillimeter wavelengths. The cosmic 
background intensity J v at z — 0, defined as the power 
received per unit frequency per unit area of detector per 
unit solid angle of sky, is given by an integral of E v (z) 
over z: 

c r a j? dt 



(21) 



We have ignored any external absorption along the line 
of sight between the sources of radiation and the ob- 
server, because all ionizing photons are assumed to be 
absorbed locally within the galaxies in which they were 
produced, and because relatively few non-ionizing pho- 
tons are absorbed by the dust in intervening galaxies. 
We have checked that, when external absorption by dust 
is included in our models, only 3% of the observed V- 
band radiation is absorbed along the line of sight from 
z = 5 to z = 0. 

2.3.2. Stellar Emissivity 
The intrinsic stellar emissivity (before absorption) 
can be expressed in terms of the comoving rate of star 
formation fL: 



(22) 



p c / dt'F„(t-t')n*(t'), 
Jo 

where F v (At) is the stellar population spectrum, defined 
as the power radiated per unit frequency per unit initial 
mass by a generation of stars with an age At. The total 
rate of star formation f2* differs from the net rate Cl s 
introduced in § 2.1 because some of the stellar material 
is ejected back to the interstellar medium by evolving 
and dying stars; thus 

(l-R)- 1 ^, (23) 



^2 j. 



where R is the IMF-averaged returned fraction. The 
observed stellar emissivity (after absorption), denoted 
by E VOl is then given by 



(1 - A V )E SV = (1 - A v )C v tl s , (24) 



where 



(25) 



,{t)J 'l-fl(f) 

is a conversion factor between the intrinsic stellar emis- 
sivity and star formation rate, which depends in general 
on the past history of (l s . At ultraviolet wavelengths, 
however, C v remains roughly constant in time, nearly 
independent of the past history of Cl s , and hence E vo (z) 
scales directly with (l s {z) in the absence of dust absorp- 
tion. 

We adopt the latest version of the Bruzual-Charlot 
population synthesis code to compute F v and R (Bruzual 
& Chariot 1998). The stellar IMF is assumed to be a 
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power law, (p(m) cx mT^ 1+x \ with x = 1.35 (a Salpeter 
IMF) and lower and upper mass cutoffs at 0.1 M Q and 
100 Mq. We further assume that all hydrogen ionizing 
photons are absorbed locally in the interstellar medium, 
where 68% of them are converted to Lya photons (the 
fraction appropriate for case B recombination in gas 
at 10 4 K). The Lya photons are assumed to be ab- 
sorbed by dust (as a result of resonant scattering by 
HI), and the remaining energy is radiated uniformly in 
wavelength between 3000 and 7000 A (a range that in- 
cludes most of the relevant emission lines). This treat- 
ment of ionizing radiation is consistent with ultravio- 
let observations of local starburst galaxies (Leitherer et 
al. 1995). In our calculations, templates of F v and R 
are pre-computed for metallicities Z/Z Q = 2.5, 1.0, 0.4, 
0.2, and 0.02. The metallicity dependence of the spec- 
tral evolution is included self-consistcntly by interpolat- 
ing the templates according to the solution for Z in our 
models of cosmic chemical evolution. Examples of these 
templates are shown in Figure 5. We will later discuss 
the effects of changes in IMF on our results (§ 4.2). 
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FIG. 5. — Spectra of stellar populations with different 
ages and metallicities. These are from the latest version 
of the Bruzual-Charlot (1998) population synthesis code 
with a Salpeter IMF truncated at 0.1 and 100 M . 

2.3.3. Absorption of Starlight by Dust 

We now consider the mean fraction A v of photons 
absorbed by the interstellar dust in galaxies. This ab- 
sorption is likely to occur on various scales of the dust 
distribution, ranging from small clumps of dust shroud- 
ing or in the vicinity of individual stars to large clouds 
in the interstellar medium. To account for all scales, it 



is appropriate first to consider the absorption along a 
ray passing through a galaxy and then to average over 
all galaxies and all impact parameters and inclination 
angles to derive the mean absorption. Along such a ray 
with fixed absorption optical depth t„, the emerging in- 
tensity is given by 



dT' v S v {T' v ) exp[-(r v - t' v % 



(26) 



where S v is the source function along the ray. In the 
absence of absorption, the corresponding intensity would 
be . Tv 

I° v {t v ) = / drX(^). (27) 



We have ignored the influence of scattering on absorp- 
tion so that S v does not contain a contribution from 
Iv To perform the average over galaxies, we define 
g(r v , z)dr v dX as the number of galaxies along a random 
ray with optical depths between r„ and r v + dr v and red- 
shift paths between X and X + dX when all galaxies and 
all positions and directions within them are considered. 
The redshift path, defined by dX = (1 + z)(l + 2q z)- 1 ^ 
for A = 0, is included in the definition of g{r v , z) so that 
an integral of T v g{j v ,z) over all t v gives the comoving 
density of dust Qd (see eq. 36). With this definition, we 
can express the mean fraction of photons absorbed by 
dust in the form 



A v {z) 



Jo° dT u g(T u ,z)[I^(T u ) - I v {t v )\ 
dT v g{T v ,z)I®{T v ) 



(28) 



The distribution function g(r V7 z) could in principle 
be modeled by specifying in detail the three-dimensional 
dust distribution in galaxies of different types, sizes, 
etc. This would, however, require a large number of as- 
sumptions. Instead, we appeal to damped Lya surveys 
to specify g{r v ,z) as a function of both r„ and z. In 
the limit of infinitely many lines of sight, these surveys 
sample correctly all impact parameters and inclination 
angles of galaxies at different redshifts. However, the 
lines of sight in quasar absorption-line surveys have fi- 
nite resolution, due to the "beam size" of quasar emis- 
sion regions. Although the exact sizes of quasar emis- 
sion regions are uncertain, the optical and X-ray con- 
tinuum regions are often assumed to be much smaller 
than 0.1 — 1 pc, which is roughly the size of the broad- 
line emission region. Thus, even in the limit of infinitely 
many lines of sight, the resulting distribution function 
from damped Lya surveys would not fully sample any 
small-scale dumpiness of the dust distribution, probably 
on scales much less than 1 pc but certainly on scales of 
clumps in the vicinity of individual stars (including cir- 
cumstellar dust). This distribution, denoted by h(r Vl z), 
has the same definition as g(r v , z), except that the for- 
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mer is an unsmoothed version of the latter. 

The relation between g and h can be written as 



r 

Jo 



d,T v p(Tl\T v )h(T v ,z), 



(29) 



where pir'J^^dT^ is the probability of finding dust 
clumps with optical depths between r' v and t' v + (Lt' v 
within a small but finite beam given that the beam has 
a mean optical depth t v . Clearly, this conditional prob- 
ability must satisfy the following relations 



/>oo 

/ dT' v p(T' v \T v ) = 1, 
JO 

dT 'v T ' v P{ T lK) = T v 



(30) 



(31) 



If the small-scale dumpiness of the dust distribution 
within galaxies were ignored, i.e., p(t' v \t v ) = <5(t^, — t u ), 
we would have g(r' v ,z) = h(r' v ^z). Thus, the probabil- 
ity distribution p introduces the additional small-scale 
dumpiness of the dust distribution specified by g but 
not sampled by h. 

We relate the distribution h of optical depths to the 
distribution / of HI column densities in the damped Lya 
systems (Fall et al. 1996). The absorption optical depth 
t v and HI column density N are related by 



mnNK a (v), 



(32) 



where n a is the mass absorption coefficient of grains and 
k m is the dust-to-HI mass ratio (eq. 14). Again, we as- 
sume that k m is only a function of rcdshift (see § 2.2.2 
and § 2.2.4), so that the distributions of t„ and N are 
simply related by h{r v , z)dT v = f(N, z)dN. From equa- 
tion (16), we then obtain 

h(T V) z) = (f t /T t )(T v /r t )-i exp(-T„/T t ), (33) 

r t (u,z) = [3cH Q /8irGf t T(2 - j)}n d (z)K a (u), (34) 

where r t is the characteristic optical depth of the ab- 
sorbers. Given the comoving density of dust Sid, the 
absorption opacity of grains K a , and the parameters f t 
(closely related to /*) and 7 discussed in § 2.2.4, the 
distribution function h is completely specified by equa- 
tions (33) and (34). It is known that the distribution 
of gas in the interstellar medium is not a perfect tracer 
of the dust distribution, especially on small scales. In 
this case, the distribution function h, when related to / 
by a constant k m , would also undersample some of the 
small-scale dumpiness of the dust in galaxies. 

The probability distribution p specifies the small-scale 
dumpiness of dust in the interstellar medium. The dust 
could be in a variety of forms (e.g., shells, spheres, and 
other complicated shapes) and in a variety of places (e.g., 
surrounding, in the vicinity of, and between stars). In- 



stead of specifying the geometry and location of such 
clumps relative to the sources of radiation, we adopt a 
statistical approach to derive p. Since we are interested 
in the mean quantity A v in galaxies, we assume as an 
idealization that all clumps at a given redshift have the 
same optical depth r c (but not necessarily the same den- 
sities, shapes, and sizes), and that such clumps are ran- 
domly distributed within the interstellar medium. The 
probability of intercepting such clumps in a small but 
finite beam passing through galaxies with a given mean 
optical depth t„ is then given by the usual Poisson dis- 
tribution: 

oc 

pK\t v ) = V ( Tl/ /T c ) n eM-r v /T c )6(T' v -nT c ). (35) 

This is similar to the clump model developed by Nata & 
Panagia (1984). It is easy to check that this probability 
distribution satisfies equations (30)-(31) and that the 
small-scale dumpiness of dust introduced by p does not 
affect the total amount of dust, i.e., 

8ttG 



n d (z) = 



3cH K a 

3cH n a 



dT v T v g{r v , z) 



dT v T v \l{T v , z) 



(36) 



Thus, regardless of how clumpy the dust is, the mean 
comoving density of dust fid computed from h is the 
same as the one computed from g. This is important 
because the mean correction factor Q derived in § 2.2.2 
is also the same regardless of whether h or g is used. 

To derive the mean fraction of photons absorbed by 
dust A v , we further assume that the source function S v 
is constant. Combining equations (28), (29), and (33)- 
(35), we obtain 



A„ = l 



1 



(7 - l)r t 



1 



— (1-e 

Tr 



(37) 



For t c — and 7=1, this reduces to the formula derived 
by Fall et al. (1996). In our new formula, the absorption 
of starlight by dust depends on the comoving density 
of dust through r t and the small-scale dumpiness of 
the dust through t c . The dependence of A„ on r t arises 
from the absorption of photons traveling through the 
whole interstellar medium when small-scale dust clumps 
are unresolved, while the dependence on t c includes the 
additional absorption by such clumps. In the limit of 
Tt S> t c , the absorption is dominated by the dust dis- 
tribution on scales resolved by quasar absorption-line 
surveys, whereas, in the limit of r t <C r c , the absorption 
is dominated by the small-scale dumpiness of the dust. 
For fixed values of r t and t c , A v depends only mildly on 
the assumption of a constant source function. For exam- 
ple, for fixed r t = 0.1 (1.0) and t c = 0, A v increases by 
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1.4% (8%) if the sources are all in the middle of the dust 
(e.g., stars embedded in clumps) and decreases by 2.6% 
(17%) if the dust is in the middle of the sources (e.g., 
clumps between stars), both relative to a constant S v . 
Therefore, the assumption of a constant source function 
is made here for simplicity and should have only a minor 
effect on A v . 

The parameter r c is by definition an absorption op- 
tical depth, and hence its dependence on frequency is 
through the absorption coefficient n a . It may also have 
some dependence on redshift as a result of the chemi- 
cal enrichment and other processes within galaxies. To 
allow for a range of possibilities, we parameterize this 
evolution by 



t c (u,z) = T c y(0) 



(") 



K a (V) 



O d (0) 



(38) 



where r c y(0) is the visual optical depth of dust clumps 
at z = and e is the power-law index that links r c to Q d . 
For e = 1, t c would evolve with redshift the same as £l d . 
This might be the case if the metals ejected by stars 
(including those locked into dust grains) were instan- 
taneously mixed in the interstellar medium, so that r c 
tracked the global chemical enrichment in galaxies. For 
e = 0, t c would remain constant in time. This might be 
the case if the timescale for mixing were longer than the 
Hubble time, so that most of the metals remained in the 
vicinity of stars with little or no evolution. Given these 
arguments, one might expect that a plausible value of e 
is somewhere between the cases of instantaneous mixing 
(e = 1) and the incomplete mixing (e = 0). We treat 
both t c v and e as adjustable parameters in our models. 

2.3.4- Dust Emissivity 

The stellar radiation absorbed by dust is reradiated 
thermally in the infrared. This emission covers a wide 
range of wavelengths, reflecting a wide range of grain 
temperatures. The far-infrared emission is dominated 
by large amounts of cold dust heated by a smooth field 
of radiation from many stars. The mid-infrared emis- 
sion is likely dominated by small amounts of warm and 
hot dust near stars, including polycyclic aromatic hy- 
drocarbons or very small grains producing some non- 
equilibrium emission. We write the dust emissivity as a 
sum of thermal emission at different temperatures: 

/•OO 

E dv = Air Pc n d K a {v) / dTB v {T)ri{T) 
Jo 

= lTT Pc £l d K a {u){B v ), (39) 

where -q{T)dT is the fraction of grains with temperatures 
between T and T + dT and (B v ) denotes the average 
of blackbodies B V {T) over the temperature distribution 
rj(T). If all grains had the same temperature T d , i.e., 



r](T) = 6(T — T d ), (By) would reduce to the spectrum 
of a single blackbody. The temperature distribution 77 
introduces the warm and hot components of the inter- 
stellar dust to account for the mid-infrared emission. 

The exact form of the dust temperature distribution 
is unknown. To specify it fully would require some 
detailed assumptions about variations of the radiation 
field, grain spatial density, and optical properties (in- 
cluding the size distribution and the grain composition) 
throughout galaxies. Since we are more interested in the 
far-infrared background, where accurate measurements 
exist, and less interested in the mid-infrared background, 
where only observational limits exist, we model the dust 
temperature distribution simply by a power law 



V(T) = 



1 IT 



6{T-T C ). 



(40) 



Here, the proportionality is fixed by the requirement 
that 77 must be normalized to unity. We have included a 
truncation factor 9(T — T c ) to exclude T < T c , because 
all grains, even those far from stars, will be heated to 
finite temperatures by the diffuse radiation field of the 
whole stellar population (or by the cosmic microwave 
background radiation if this exceeds the diffuse starlight 
in galaxies). Since the total emission must be finite, the 
dust temperature index in equation (40) must satisfy 
a > 5 + n for K a (v) oc v n . Indeed, the observed mean 
infrared colors of nearby galaxies can be fitted with a 
value of a that satisfies this requirement (§ 2.3.5). For 
simplicity, we have not imposed any upper cutoff tem- 
perature on 77, since this has little effect on our results 
when a is large. Such a steep slope is the consequence 
of a very small fraction of mass in warm and hot dust in 
the interstellar medium. 

The lower cutoff temperature T c in equation (40) is 
closely related to a mean temperature of dust T d , defined 
here by the first moment of T over r)(T): 

T d = dTTi](T) = (^—^j T c . (41) 

Given the dust temperature index a, the averaged black- 
body spectrum (B v ) is then a unique function of the 
mean dust temperature T d . We determine the redshift 
dependence of this mean temperature by the usual con- 
dition of energy balance 



4np c n d / dvK a {v)[{B v ){T d ) - B„(T C mb 
Jo 



)] 



dvA v E s 



(42) 



where T CMB = 2.728(1 + z) K is the temperature of 
the cosmic microwave background radiation (Fixsen et 
al. 1996). The inclusion of the cosmic microwave back- 
ground has little effect on T d at z <J 7 in our models. 
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We have ignored other sources of radiation (primarily 
active galactic nuclei) in the heating of interstellar dust. 
The blue emissivity of quasars, estimated from the ob- 
served luminosity function (Pei 1995), is less than a few 
percent of the blue emissivity of galaxies at all redshifts 
in our models. The far-infrared emissivity of quasars, 
however, is not yet known. At long wavelengths [hv <C 
kTd), the multi-temperature spectrum {B v )(Td) reduces 
to the Rayleigh- Jeans part of a single blackbody with the 
same Td- At short wavelengths (hv 3> kTd), the mean 
spectrum exceeds the Wien tail of a single blackbody at 
the mean temperature. 

2. 3. 5. Inputs from Galaxy Surveys 

The formulae developed above allow us to input the 
observed emissivity E vo from optical imaging and red- 
shift surveys into our models of cosmic chemical evolu- 
tion. The observable E vo specifies the comoving rate of 
star formation f2 s through the relation 



1 



l-A v (z) 



C v {z) 



(43) 



where C v converts between the intrinsic stellar emissiv- 
ity and star formation rate and A v corrects for the stel- 
lar photons absorbed by dust in the interstellar medium 
(eq. 24). Because A v is coupled to fid (eqs. 34 and 
37) and C v is coupled to fl s (eq. 25), solutions to our 
models of cosmic chemical evolution must be derived 
iteratively (see § 3). In the following, we adopt the ob- 
served ultraviolet emissivity E vo {z) from optical surveys 
of galaxies and estimate the dust temperature index a 
from observations of IRAS galaxies in the local universe. 
This then leaves two unknown parameters in our models: 
t c v mainly affects the overall amplitude of the absorp- 
tion by dust and hence the amplitude of the far-infrared 
background intensity, whereas e mainly affects the evo- 
lution of the dust opacity and hence the spectral shape 
of the background intensity. Thus, T c y and e can be 
constrained by measurements of the cosmic far-infrared 
background. 

Estimates of the observed ultraviolet emissivity E uo (z) 
are now available over the redshift interval 0.2 <J z <^ 4.5, 
from the Canada- France Redshift Survey (CFRS, Lilly 
et al. 1995) and the Hubble Deep Field (HDF, Williams 
et al. 1996). Table 2 lists the estimates of E vo {z) at 
several redshifts for an assumed cosmology with A = 0, 
go = 1/2, and H = 50 km s^ 1 Mpc -1 . The entries at 
z = 0.35 - 0.875 were derived by Lilly et al. (1996) from 
the rest-frame 2800 A luminosities of about 600 galaxies 
in the CFRS sample with spectroscopically confirmed 
redshifts over < z < 1. The entries at z = 1.25 — 1.75 
were derived by Connolly et al. (1997) from the rest- 
frame 2800 A luminosities of some 200 galaxies selected 



in the HDF images with photometrically estimated red- 
shifts over 1 < z < 2. The two entries at z = 2.75 
and z = 4.0 were derived by Madau et al. (1996, 1998): 
the former is based on 69 ultraviolet-dropout galaxies 
in the HDF images with color-estimated redshifts of 
2.0 <; z <J 3.5, and the latter is based on 14 blue dropout 
galaxies with color-estimated redshifts of 3.5 <J z <J 4.5. 
The ultraviolet emissivity at z — is not directly ob- 
served. Instead, we use the mean B-band luminosity 
density of local galaxies from Ellis et al. (1996), listed as 
the first entry of E vo (z) in Table 2. The values of E vo 
derived by these authors all include modest extrapola- 
tions beyond the observed range of luminosities using a 
Schechter function with a « —1.3. The quoted statis- 
tical uncertainties in E vo are 20 — 30% at z <. 1 and 
40 - 60% at z ;> 1. 

We determine the temperature index a of the dust 
by fits to the observed infrared emissivity of the local 
universe from the IRAS survey. The results are shown 
in Figure 6, where £^(0) is normalized at Ai = 100 /im, 
so that the ratio of Edu(0)/Ed Vl (0) depends mainly on 
a. The data points were derived by Soifer & Ncugc- 
bauer (1991) from complete flux-limited samples (z < 
0.08 and (z) = 0.006) at A = 12, 25, 60, and 100 /xm. 
The plotted errors reflect statistical uncertainties in the 
corresponding luminosity functions tabulated by Soifer 
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FIG. 6. — Normalized infrared emissivity at z = 0. 
The data points were derived by Soifer & Neugebauer 
(1991) from IRAS observations of nearby galaxies at 
A = 12, 25, 60, and 100 /im. The solid curves represent 
the multi-temperature models of dust emission with the 
indicated values of the temperature index a, while the 
dashed curve represents a single blackbody. 
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& Neugebauer (1991). The solid curves show the multi- 
temperature models with the three indicated values of a, 
while the dashed curve shows a single blackbody (with 
a = oo). All curves are computed with K a oc v 2 , an ap- 
proximation at A > 10 (/m for LMC-type dust (§ 2.2.5), 
and Td(0) = 16.5 K, a typical mean temperature of 
dust at z = obtained in our models (§ 3.3), but these 
curves are insensitive to the exact value of T^, varying by 
roughly 30% at 12 /xm for T d w 10-25 K. In the context 
of this work, namely, using the cosmic infrared back- 
ground to correct for any missing starlight in the optical 
data, it is important to have a good account of the total 
infrared energy of the dust emission. A single blackbody, 
if adopted, would underestimate the mid-infrared emis- 
sion, missing roughly 50% of the total infrared energy. 
We adopt the multi-temperature spectrum with a — 8.1 
in all of our calculations. This still underestimates the 
emission near 12 /xm. The observed excess is believed to 
be contributed by very small grains (including polycyclic 
aromatic hydrocarbons), which can be excited by sin- 
gle photon hits and produce non-thermal emission with 
complex emission features (Desert et al. 1990). Never- 



theless, the mid-infrared excess near 12 /xm contributes 
less than 10% of the total infrared energy and hence 
should have little effect on our fits to the cosmic infrared 
background. 

3. RESULTS 

Our models of cosmic chemical evolution are specified 
by equations (4)- (6). These equations link the history 
of metal enrichment to the histories of star formation 
and gas consumption in galaxies and the baryon flow 
between the interstellar and intergalactic media. We 
adopt two observational inputs to solve the models: (1) 
the observed comoving density of HI, Ohio, from damped 
Lya surveys, and (2) the observed ultraviolet emissivity 
of galaxies, E VOl from optical imaging and redshift sur- 
veys. The first is related to the comoving density of gas, 
Sl g , (eq. 19), including a correction factor, Q, for the ab- 
sorbers missing from existing damped Lya surveys that 
is self-consistently coupled to the comoving density of 
dust (eq. 17). The second is related to the comoving 
rate of star formation, Q s , (eq. 43), which depends on 



TABLE 1 
Model Inputs 



Parameter 


Value 


Source 


Present mean metallicity 


Z(0) = Z Q = 0.02 


solar value 


Mean dust-to-metals ratio 


d m = 0.45 ±0.1 


fit - Fig. 1 


Slope of bright quasar counts 


= 2.0 ±0.1 


fit - Fig. 2 


Index of gamma distribution 


7 = 1.2 ±0.1 


fit - Fig. 3 


Norm of gamma distribution 


/* = 0.03 ± 0.01 


fit - Fig. 3 


Index of temperature distribution 


a = 8.1 ±0.3 


fit - Fig. 6 


Opacities of grains K e & n a 


LMC-type dust 


model - Fig. 4 


Population synthesis F v & R 


Salpeter IMF, 0.1 - 100 M© 


model - Fig. 5 


Observed density of HI $lw (z) 


Data - absorption surveys 


listed in Tab. 2 


Observed UV emissivity E uo (z) 


Data - optical surveys 


listed in Tab. 2 




f (1.13,1.0) 


minimum J„ fit 


Small-scale clumpy ISM model 


(r cV ,e) = \ (1.06,0.4) 


best J„ fit 




L (1.13,0.1) 


maximum J„ fit 



NOTE. — The present mean metallicity Z(0) is taken to be the solar value (§ 2.1). The dust-to-metals mass ratio d m 
is the mean value in nearby galaxies and damped Lya systems at 0.7 < z < 2.8 (§ 2.1). The bright-end slope (3 of 
the luminosity function of quasars is fitted from the observed counts at B < 20 (§ 2.2.3). The gamma distribution 
parameters 7 and /» for the distribution of HI column densities are fitted from surveys of damped Lya systems 
(§ 2.2.3). The dust temperature index a is estimated from the mean infrared colors of local galaxies (§ 2.3.5). The 
extinction and absorption opacities of dust grains, n e and K a , are derived from the graphite-silicates model that fits 
the observed LMC extinction curve (§ 2.2.5). The population synthesis spectrum F v and returned fraction R are 
computed from the latest version of the Bruzual-Charlot models with a Salpeter IMF truncated at 0.1 and 100 M 
(§ 2.3.2). The present visual optical depth T c y of small-scale dust clumps and the evolution index e of the clumps 
relative to the global chemical enrichment history of galaxies are determined from fits to the DIRBE and FIRAS 
measurements of the cosmic infrared background (§ 3.4). 



14 



a conversion factor, C v , between the intrinsic emissivity 
and star formation rate (with the aid of stellar popula- 
tion synthesis models in eq. 25) and a correction factor, 
A v , for stellar radiation absorbed by dust within galaxies 
that is again self-consistently coupled to the comoving 
density of dust (eq. 37). The background intensity, J„, 
at the present epoch is given by equation (21), including 
the radiation from both stars (eq. 24) and dust (eq. 39). 

Table 1 lists all the input parameters of our models. 
The only adjustable parameters are the present visual 
optical depth T c y of small-scale clumps of dust and the 
index e that specifies the evolution of these clumps rel- 
ative to the mean comoving density of dust (§ 2.3.3). 
Since the obscuration, absorption, and emission by dust 
are interrelated, all the equations in our models are cou- 
pled. Iterative solutions are obtained as follows. Given 
t c v and e, the initial set of Q g and £l s (and hence all 
other quantities as functions of rcdshift in our models) 
is computed with Q = 1, A v = 0, and C v from a con- 
stant tt s . The values of Q, A v , and C v are then updated 
from this initial solution so that the next solution can 
be computed. This iterative process is repeated until 
fl g and Cl s (and hence all the other quantities) converge 
at all redshifts. The stellar and dust emissivities, E sv 



and Ed„, as functions of z, and the background inten- 
sity, J,,, at z = are then computed from the converged 
solutions. Fits of the far- infrared part of J v to the back- 
ground measurements then give the values of T c y and e. 
These parameters are not strongly correlated; r c y de- 
termines mostly the amplitude of the far-infrared back- 
ground, while e controls mostly the spectral shape of the 
background. In addition to the best-fit solution, we also 
present the 95% confidence minimum and maximum so- 
lutions to bracket the observational uncertainties in the 
infrared background measurements (described in detail 
in § 3.4). The corresponding values of r c y and e arc 
listed in Table 1. 

3.1. Interstellar Gas and Star Formation 

The comoving density of interstellar gas Q g and the 
comoving rate of star formation (l s are the two key quan- 
tities from which all other quantities in our models can 
be calculated. Table 2 lists selected solutions (f2 g , Z, 
Cl s , A v ), together with the input data flmo and E vo . 
We have adopted discrete input data, because flmo(z) at 
different z is obtained from samples of quasars selected 
at different effective wavelengths A e , while E vo {z) at dif- 



TABLE 2 
Selected Solutions 



z 




A e 




Q g (z) 






Z(z)/Z Q 


0.00 


(0.25) 


5500 


0.50 ±0.12 




(1.000) 


0.64 


0.54 ±0.31 


5500 


2.28 


(1.90, 


2.45) 


0.511 


(0.473, 0.524) 


1.89 


1.58 ±0.92 


5500 


6.58 


(4.04, 


8.43) 


0.102 


(0.076, 0.111) 


2.40 


2.15 ±0.70 


5500 


6.88 


(4.23, 


9.68) 


0.060 


(0.036, 0.070) 


3.17 


2.34 ± 1.14 


5500 


4.66 


(3.65, 


6.65) 


0.026 


(0.013, 0.040) 


4.01 


1.44 ±0.58 


6500 


2.15 


(1.99, 


2.63) 


0.016 


(0.008, 0.035) 




log E vo {z) 


A 




ti a {z) 






A v {z) 


0.000 


19.47 ±0.10 


4400 


0.60 


(0.66, 


0.63) 


0.468 


(0.486, 0.487) 


0.350 


18.89 ±0.07 


2800 


1.97 


(2.63, 


1.86) 


0.645 


(0.718, 0.627) 


0.625 


19.21 ±0.08 


2800 


4.88 


(6.52, 


4.23) 


0.688 


(0.758, 0.650) 


0.875 


19.53 ±0.15 


2800 


11.2 


(15.0, 


9.66) 


0.691 


(0.760, 0.651) 


1.250 


19.69 ±0.15 


2800 


17.0 


(19.9, 


15.3) 


0.671 


(0.704, 0.642) 


1.750 


19.59 ±0.15 


2800 


11.1 


(9.13, 


10.6) 


0.630 


(0.507, 0.635) 


2.750 


19.42 ±0.15 


1500 


14.7 


(4.72, 


20.5) 


0.823 


(0.472, 0.874) 


4.000 


19.02 ±0.20 


1500 


2.97 


(1.21, 


7.29) 


0.637 


(0.118, 0.848) 



NOTE. — Input quantities are listed in the form observation ± la uncertainty Output quantities are listed in the 
form best (minimum, maximum) solutions allowed by the cosmic far-infrared background. All entries pertain to 
A = 0, qo = 1/2, and Ho — 50 km s _1 Mpc -1 ; Qmo and Q g in units of 10~ 3 , A e and A in units of A, E vo in units of 
W Hz -1 Mpc -3 , (l s in units of 10~ 4 Gyr -1 . The relation between the total and net star formation rates is Q*(z) = 
(1 - Ry-n^z), with R = 0.33, 0.32, 0.31, 0.30, and 0.28 at z = 0, 1, 2, 3, and 4, respectively. 
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FIG. 7. — Comoving density of interstellar gas (upper panel) and comoving rate of star formation (lower panel). 
The filled circles are the pointwise solutions, connected by the solid, short-dash, and long-dash curves to distinguish 
the best, minimum, and maximum solutions, respectively. In the upper panel, the open circles with error bars, 
connected by the dotted curve, are the observed comoving density of neutral gas in damped Lya systems. In the 
lower panel, the open circles with error bars, connected by the dotted curve, are the observed ultraviolet emissivities 
when expressed as star formation rates (without correction for absorption by dust). 



ferent z is measured at different rest-frame wavelengths 
A. Accordingly, Q(z) is evaluated at the same redshifts 
as Ohio (z), while A v (z) and C v {z) (averaged through a 
±10% wavelength box filter) are evaluated at the same 
redshifts as E uo (z). The discrete values of Q(z), A v (z), 
and C v (z), along with Clmo(z) and E uo {z), are then suf- 
ficient to define the pointwise solutions for Sl g (z) and 



Sl 3 (z). To derive other functions of redshift, we inter- 
polate in ^lg(z) and Sl a (z) using the form logfi(z) oc z. 
With other interpolation schemes, such as logf2(z) oc 
log(l + z), the results would differ typically by a few 
percent, reaching only 10% at z « 4. 

Figure 7 shows the comoving density of interstellar 
gas and the comoving rate of star formation as functions 
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of redshift. The filled circles are the model solutions for 
£l g {z) and fl s (z), connected by the solid, short-dash, 
and long-dash curves to distinguish the best, minimum, 
and maximum (95% confidence) solutions, respectively. 
In the upper panel, the open circles with error bars are 
the observed comoving density of neutral gas in damped 
Lya systems, /if^Hio- The obscuration correction to Q g 
(from the open to filled circles) isQ = 1.5,2.5 — 3.7, and 
1.3 — 2.5 at z = 0, 1, and 3, where the ranges include the 
three solutions. This correction is large at z w 1, when 
the comoving density of dust is near its highest value. 
The peak in tt g occurs at z k 2.4, just before gas is 
consumed faster by star formation than it is replenished 
by inflow. In the lower panel, the open circles with error 
bars are the observed (i.e., optically inferred) star for- 
mation rates derived in the absence of dust absorption, 
E vo /C v . The absorption correction to Q s (from the open 
to filled circles) is (1 - A,,)- 1 = 2.4 - 2.5, 2.8 - 3.9, and 
1.7 — 7.8 at z = 0, 1, and 3. All three solutions for Cl s (z) 
show a rapid decrease with respect to time at z <J 1 and 
an increase with time at z ;> 3. The best and maximum 
solutions remain roughly flat over 1 <; z <; 3, while the 
minimum solution still shows an increase with time over 
this redshift range. 

The star formation rates in our models differ among 
the three solutions by less than 30% at z <; 2 but by a 
factor of three at z ;> 3. The similarity at low redshifts 
is a consequence of the constraint imposed by the far- 
infrared background. In our models, the mean comoving 
density of dust is small at high redshifts, and hence sub- 
stantial absorption occurs only if the dust is extremely 
clumpy, as would be the case for heavily embedded stars. 
Indeed, we find that the minimum solution corresponds 
to an instantaneous mixing of dust in the interstellar 
medium (i.e., small-scale clumps follow the evolution of 
the whole interstellar medium), whereas the maximum 
solution corresponds roughly to an incomplete mixing 
of dust (i.e., small-scale clumps are nearly independent 
of the evolution of the rest of the interstellar medium) , 
with the best solution representing an intermediate case. 
Therefore, the different solutions derived here reflect dif- 
ferent evolution of the small-scale inhomogeneity of dust 
in galaxies. 

3.2. Chemical Enrichment of the Interstellar Medium 

Figure 8 shows the mean abundance, Z, and comov- 
ing density, Q m , of heavy elements in the interstellar 
medium as functions of redshift in our models. Again, 
the solid, short-dash, and long-dash curves represent the 
best, minimum, and maximum solutions, respectively. 
The data points at z > 1.5 were derived by Pettini 
et al. (1997b) from the observed abundances of Zn in 



damped Lya systems, while the data point at z = 0.77 
was derived by Boisse et al. (1998) by combining their 
own observations with those of Pettini et al. (1997b). 
These are the HI column density weighted mean metal- 
licity and hence are the correct measure of metallicity to 
be compared with our models. The uncertainty in the 
data is still large and no attempt is made here to cor- 
rect for a preferential selection of the low-metallicity ab- 
sorbers that would bias the observed mean. We empha- 
size, however, that this bias could potentially be large at 
z ~ 1, where the comoving density of dust is high. Nev- 
ertheless, all three solutions for the mean metallicity in 
our models are consistent within the uncertainties with 
the observations, especially if the data point at z = 0.77 
is regarded as a lower limit. The present mean metal- 
licity Z(0), taken here to be the solar value, implies an 
effective yield of y = 0.45Z Q , nearly the same in the 
three solutions. The peak in f2 m at z « 0.8 is due to a 
combination of the rise in Z and the decline in fl g toward 
z = 0. The decrease in Q m toward z = is caused partly 
by the fact that star formation starts to consume more 
metals than it contributes to the interstellar medium (at 
z <; 0.8, when Z exceeds y), and partly by some outflow 
of the interstellar metals to the intergalactic medium 
(see § 4.4). 




z 



FIG. 8. — Mean metallicity of interstellar gas in galax- 
ies in units of the solar value (upper panel) and comov- 
ing density of heavy elements in galaxies (lower panel) . 
The solid, short-dash, and long-dash curves are the best, 
minimum, and maximum solutions in our models. In the 
upper panel, the data points at z > 1.5 are from Pettini 
et al. (1997b), while the data point at z = 0.77 is from 
Boisse et al. (1998). 
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FIG. 9. — Comoving emissivity E u (z) times frequency v as a function of redshift. The solid, short-dash, and long- 
dash curves are the best, minimum, and maximum solutions in our models. The data points are from Madau et 
al. (1998) for 0.15 ^m; Lilly et al. (1996) for 0.28, 0.44, and 1.0 fim at z = 0.25-1.0; Connolly et al. (1997) for 0.28 
and 0.44 /um at z — 1 — 2; Ellis et al. (1996) for 0.44 /im; Gardner et al. (1997) for 2.2 /im; and Soifer & Neugebauer 
(1991) for 25, 60, and 100 /im. All wavelengths are in the rest frame at the indicated redshifts. 



The histories of metal enrichment in the three solu- 
tions presented here differ mainly at high redshifts and 
are caused by the different histories of star formation. 
This indicates that the observed mean metallicity in the 
high-rcdshift damped Lya systems should be a useful 
constraint on star formation at high redshifts. In par- 
ticular, if the global rates of star formation were much 
higher at z ;> 3 than in the maximum solution presented 



here, the interstellar metallicity would be larger than 
that observed in the damped Lya systems, unless a sig- 
nificant outflow had occurred before z ~ 3. 

3.3. Evolution of Stellar and Dust Emissivities 

Figure 9 shows the emissivity E v {z) times frequency 
v as a function of redshift at several rest-frame wave- 
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lengths. The data are taken from Lilly et al. (1996) for 
0.28, 0.44, and 1.0 /im at z = 0.35, 0.625, and 0.875, 
Madau ct al. (1998) for 0.15 at z = 2.75 and 4.0, 
Connolly et al. (1997) for 0.28 and 0.44 jimatz = 1.25 
and 1.75, Ellis et al. (1996) for 0.44 patz = 0, Gard- 
ner et al. (1997) for 2.2 at z = 0, and Soifcr & 
Neugebaucr (1991) for 25, 60, and 100 fim at z = 0. 
The curves are the three solutions for the cosmic emis- 
sivity in our models. The results at A = 0.15 — 2.2 fim 
represent the direct but attenuated starlight, while those 
at A = 25 — 100 /urn represent the starlight reradiated 
by dust. The data that were input to the models are 
plotted as open circles, while those that are output are 
plotted as filled circles. Evidently, our models reproduce 
remarkably well even data that were not used as input. 
In particular, the agreement at A = 1.0 — 2.2 /im, where 
the absorption by dust is relatively small, is an indica- 
tion that the Salpeter IMF is reasonable. A steeper IMF, 
such as a Scalo IMF, would produce too much 2.2 /im 
light at z = (by a factor of two). The agreement with 
the observed 25 — 100 /im cmissivity at z = indicates 
that the evolution of the small-scale dumpiness of dust 
in the models is reasonable. A substantially different 
evolution from that presented here might give an accept- 
able fit to the integrated infrared background but would 
not necessarily match the local infrared emissivity at the 
same time. 

The mean temperature of dust, plotted in the upper 
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FIG. 10. — Mean temperature of interstellar dust (up- 
per panel) and ratio of bolometric dust emissivity to 
bolomctric stellar emissivity (lower panel). The solid, 
short-dash, and long-dash curves represent the best, 
minimum, and maximum solutions, respectively. 



panel of Figure 10, varies from Td = 21 — 27 K at z = 
4 to Td — 22 — 24 K at z — 2 and then decreases to 
Td = 16.5 K at z — 0. In our models, dust grains heated 
to different temperatures (at least above a finite tem- 
perature T c w 0.9Td) contribute different amounts to 
the bolomctric dust emission: 50% by dust with T « 
(0.9 - 1.6)T d , 25% by dust with T w (1.6 - 3.0)T d , and 
25% by dust with T > 3.0T d , with little dependence 
on redshift. The ratio of the bolometric dust emission 
to the bolomctric stellar emission, defined by Ed/E s — 
J dvEdvj J dvE sv and plotted in the lower panel of Fig- 
ure 10, varies from Ed/E s = 14 - 68% at z = 4 to 
Ed/E s = 54 — 64% at z — 2 and then decreases to 
Ed/ E s — 36% at z = 0. Evidently, the dust absorbs 
and reradiates a significant fraction of the stellar en- 
ergy. Moreover, since r c is typically larger than r t , most 
of this absorption occurs in small-scale clumps, rather 
than the large-scale distribution of dust in the interstel- 
lar medium. This is consistent with the general notion 
that the young stars responsible for most of the energy 
production arc typically located in dusty star-forming 
regions. 

3.4. Cosmic Background Radiation 

Figure 11 shows the background intensity J„ at z = 
times frequency v as a function of wavelength A. The 
circles at 140 and 240 ^m are the DIRBE detections; 
the errors include both statistical and systematic uncer- 
tainties in the measurements and the removal of fore- 
ground emissions (Hauser et al. 1998). The zigzag curve 
at 150 - 2000 ^m is the weighted mean of the FIRAS 
detections, based on three independent methods to re- 
move the foregrounds (Fixsen et al. 1998). The statis- 
tical uncertainties in the FIRAS results are 40%, 27%, 
5%, 6%, 10%, and 15% at A = 160, 180, 240, 400, 700, 
and 1000 fim, respectively (Fixsen 1998, private com- 
munication). The solid curve is the weighted minimum- 
X 2 fit of the background intensity in our models to the 
far-infrared/submillimcter DIRBE and FIRAS measure- 
ments to determine T c y and e. We regard this as our 
"best" solution. The systematic uncertainties in the FI- 
RAS results are not known precisely. For this reason, 
two other acceptable fits are presented as the short-dash 
and long-dash curves. These "minimum" and "maxi- 
mum" solutions would be the 95% confidence lower and 
upper limits if the systematic uncertainties in the FIRAS 
data were taken to be the difference between the three 
independent methods in the removal of the foregrounds 
by Fixsen et al. (1998). 

Also shown in Figure 1 1 is a comparison of the back- 
grounds in our models with various measurements and 
limits at other wavelengths. The squares with errors 
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FIG. 11. — Extragalactic background intensity J„ times frequency v as a function of wavelength A. The circles are 
the DIRBE detections of the cosmic infrared background at 140 and 240 /im (Hauser et al. 1998). The zigzag curve 
is the FIRAS detections at 150 — 2000 /im (Fixsen et al. 1998). The solid, short-dash, and long-dash curves represent 
the best, minimum, and maximum fits of our models to these observations. The filled squares are the integrated 
background light from galaxy counts at 2000 A (Armand et al. 1994) and 3600 - 22000 A (Pozzetti et al. 1998). The 
open circle is a tentative detection of the extragalactic background at 3.5 fim (Dwek & Arendt 1998). The arrows 
and the stepped line indicate various observational limits on the extragalactic background radiation (see § 3.4 for 
references) . 



are the integrated light from observed sources at 2000 A 
(Armand, Milliard, Deharveng 1994) and 3600-22000 A 
(Pozzetti et al. 1998). The large error at 2000 A is due to 
the extrapolation of the counts beyond a limiting mag- 
nitude of 18.5. At optical wavelengths, the slope of the 



counts suggests that very faint uncounted galaxies con- 
tribute little to the background. Thus, we regard the in- 
tegrated light from the observed sources as tantamount 
to a measurement of rather than a lower limit on the 
background. The open circle is a tentative detection of 
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the extragalactic background at A = 3.5 fim from the 
DIRBE (2.2 and 3.5 /xm) maps and an assumed value 
of the background at A = 2.2 ^im by Dwek & Arendt 
(1998). Various (downward and upward) arrows are ob- 
servational (upper and lower) limits: the DIRBE 2a up- 
per limits at 1.25 — 100 /xm (Hauser et al. 1998), the high 
galactic latitude upper limits at 1000 A (Holberg 1986) 
and 1700 A (Bowyer 1991), the sky surface photometry 
upper limits at 4400 A (Toller 1983) and 5100 A (Dube, 
Wickes, & Wilkinson 1979), the lower limits in the no- 
evolution case from IRAS galaxy counts at 25, 60, and 
100 /im (Hacking & Soifer 1991), the upper limits at 
10 — 100 fim from the fluctuation analysis of the DIRBE 
maps (Kashlinsky, Mather, & Odenwald 1996), and the 
upper limits at 5 — 25 /im (Stanev & Franceschini 1998) 
and at 1 — 50 ^m (stepped line, Biller et al. 1998) from 
the lack of attenuation in the TeV gamma-ray spectrum 
of Mrk 501. Evidently, our models agree well with a wide 
variety of measurements of and limits on the background 
intensity over four decades in wavelengths. 

The total energy density radiated by stars of all ages 
and now contained in the extragalactic background (over 
912 A < A < 5 mm) is, in units of the present critical 
energy density, 

tt R =— / dv.J v = (5.1 -5.5) x HT 6 . (44) 

Of this, 54 — 59% is reprocessed by dust into the cosmic 
infrared background. For comparison, the energy den- 
sity in the cosmic microwave background radiation left 
by the Big Bang is about 20 times larger than that given 
by equation (44) (Fixsen et al. 1996). In these estimates, 
we have ignored the unknown contribution from active 
galactic nuclei. 

4. DISCUSSION 

We have combined our models with a wide variety of 
observations to study the average evolution of the entire 
population of galaxies. Before discussing some implica- 
tions of our solutions (§ 4.4 and § 4.5), we examine some 
major uncertainties that could potentially affect our re- 
sults (§ 4.1 and § 4.2) and compare them with other 
estimates (§ 4.3). 

All of the results presented here are based on a cos- 
mological model with A = 0, qo = 1/2, and H = 
50 km s _1 Mpc -1 . These parameters may no longer be 
regarded as "standard" in light of recent evidence, espe- 
cially from high-redshift supernovae (see, for example, 
Reiss et al. 1998). However, many of the observational 
input data on which our results depend were reported 
only for this particular set of cosmological parameters. 
A self-consistent exploration of the effects of different 



cosmological parameters would therefore require a re- 
analysis of the input data, particularly Ohio and E uo . 
Such a study might prove interesting but is beyond the 
scope of this paper. 

4.1. Model Uncertainties 

Our results are independent of many properties of the 
damped Lya systems, including their morphologies, pro- 
vided only that they represent the main reservoirs of gas 
for star formation. The existing evidence, although not 
yet conclusive, suggests that the damped Lya systems 
are galaxies of a wide variety of morphological types (Le 
Brun et al. 1997). The idea that the damped Lya sys- 
tems trace the the bulk of star formation is supported 
by the observation that they contain enough neutral gas 
at z * 2 — 3 to make all the stars visible today (Wolfe 
et al. 1995). Such a large amount of neutral gas is not 
present in the local universe, even including the con- 
tribution of low surface brightness and dwarf galaxies 
(Zwaan et al. 1997). Indeed, if the neutral gas in damped 
Lya systems were not converted into stars, it must have 
been expelled or ionized at low redshifts. 

The main approximation in our treatment of the in- 
terstellar gas in galaxies is the neglect of molecules. As 
discussed in § 2.2.1, this could increase f2 s by up to a 
factor of two, most likely at low redshifts. Accordingly, 
would be also be higher because it is linked to the 
metallicity and the present value Z(0) is fixed in our ap- 
proach. However, as we have already emphasized, Z(0) 
is not known precisely, and different values of this pa- 
rameter affect fid through Q g as well. We have checked 
that these uncertainties, comparable to the errors in the 
observed comoving density of HI at low redshifts, do not 
strongly affect our solutions for the global rate of star 
formation. We have also checked that the uncertainties 
in the input parameters (d m , f3, 7, /*, a) all have minor 
effects on our results. 

A natural consequence of dust in galaxies is that it 
obscures background quasars (Ostriker & Heisler 1984). 
This effect also causes incompleteness in samples of 
damped Lya systems found in the spectra of optically se- 
lected quasars (Fall & Pei 1993). If this bias were ignored 
in our models, the resulting mean interstellar metallic- 
ity at z « 2.5 would be 3 — 5 times higher than that 
observed in the damped Lya systems. As discussed at 
length in § 2.2.2, our models include self-consistent cor- 
rections for missing damped Lya systems and quasars. 
From our solutions for / and f Q and the formulae of 
Fall & Pei (1993), we estimate that the missing frac- 
tion of optically selected quasars with V <^ 20 (R <^ 
20) is relatively small: 9 - 12% (8 - 11%), 12 - 19% 
(10 - 17%), and 13 - 23% (12 - 21%) at z = 2, 3, and 
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4, respectively. An important observational constraint 
on obscuration comes from counts of radio and optically 
selected quasars on the assumption that they have the 
same intrinsic evolution. Shaver et al. (1996) found that 
the comoving densities of radio and optically selected 
quasars have similar apparent evolution to z « 4. How- 
ever, the samples at high redshift are very small and are 
easily consistent with 40 — 50% of the optically selected 
quasars being missed at z w 3 — 4 (see Fig. 2 of Shaver et 
al. 1996). The obscuration in our models is well within 
this limit. 

The determination of the global star formation rate 
from the observed ultraviolet emissivity of galaxies in- 
volves assumptions about the IMF and about dust 
(eq. 43). In our approach, all uncertainties in the IMF 
(its shape, mass cutoffs, and any redshift evolution) are 
combined in the conversion factor C v {z) between the 
intrinsic emissivity and star formation rate (§ 2.3.2). 
All uncertainties related to the dust (its opacity, evo- 
lution, and distribution) are combined in the mean frac- 
tion A v (z) of photons absorbed by dust (§ 2.3.3). Our 
solutions are based on the following assumptions: (1) an 
IMF that is constant in time; (2) an IMF with Salpetcr 
form; (3) a constant dust-to-gas ratio and a power-law 
relation between the optical depth of small-scale dumpi- 
ness and mean comoving density of dust; (4) optical 
opacity of grains as in LMC-type dust; and (5) infrared 
opacity as in LMC-type dust. We now discuss each as- 
sumption in turn. The basic conclusion is that the small- 
scale dumpiness of the dust and the optical opacity of 
grains are the major sources of uncertainty in the deter- 
mination of the global star formation history. 

(1) Some of the main inputs to our models are the ob- 
served rest-frame 1500 A and 2800 A emissivities. These 
models, with a constant IMF, also reproduce as output 
the observed rest-frame 4400 A and 1 /im emissivities at 
0.2 <; z <; 2.0 (Fig. 9). This may be an indication that 
any evolution of the IMF is relatively weak, at least for 
z <J 2. Since absorption by dust decreases with wave- 
length, the rest-frame 2.2 /im emissivity of galaxies, if 
observed over a similar range of redshifts, could provide 
a tighter constraint on this evolution than the rest-frame 
4400 A and 1 fj,m emissivities. 

(2) The Salpeter IMF is flatter (has more massive 
stars) than the Scalo IMF, the other IMF often employed 
in studies of stellar populations. The appropriateness of 
the Salpeter IMF over the Scalo IMF to reproduce the 
observed emissivities at z — has already been noted by 
Lilly et al. (1996) and Madau et al. (1998). Repeating 
all of our calculations with a Scalo IMF, we find that 
the solutions for Cl s (z) are not strongly affected, except 
at very low redshifts (z <^ 0.5). However, a Scalo IMF 
would predict too much red light, causing a factor of 



roughly 2 — 3 mismatch of our models with the observed 
2.2 /im local emissivity and with the 1 — 2.2 /j,m extra- 
galactic background. We have not fine-tuned the upper 
and lower mass cutoffs of the IMF. The lower mass cut- 
off could potentially affect the overall amplitude of the 
star formation rate, without any strong effects on the 
stellar emissivity. With a lower mass cutoff of 0.1 M Q , 
the comoving density of stars at z = in our models 
agrees well with the empirical estimates (§ 4.3). 

(3) Our models of cosmic chemical evolution repre- 
sent a self-consistent way to relate the evolution of dust 
to the history of star formation on the assumption of a 
constant dust-to-metals ratio. If this ratio were to vary 
by 50% over <; z <; 3 (Fig. 1), it would have a mi- 
nor effect on the solutions for fl s (z). We have adopted 
a simple statistical approach to characterize the small- 
scale dumpiness of dust in the interstellar medium. This 
introduces two adjustable parameters: the present value 
of the optical depth of the clumps r c and an index e that 
specifies how the clumps evolve relative to the mean co- 
moving density of dust. In our models, the small-scale 
dumpiness of dust is responsible for most of the ab- 
sorption of starlight. We have presented three solutions 
in which the small-scale dumpiness of the dust evolve 
differently. This turns out to be a major source of un- 
certainty in our solutions for the global star formation 
rate. 

(4) The optical opacity of grains varies between galax- 
ies by a factor of two, even among the Milky Way, 
LMC, and SMC (Fig. 4). To assess this effect on our 
results, we repeated our calculations with Galactic and 
SMC- type dust. The resulting infrared background was 
similar within 30%, while the ultraviolet background at 
A « 2000 A differed by a factor of two higher (lower) 
for Galactic-type (SMC-type) dust relative to the LMC- 
type dust. The emissivity E vo at rest-frame ultraviolet 
wavelengths increased (decreased) by a factor of two at 
low redshifts and decreased (increased) at high redshifts. 
The solution for (l s (z) was similar at z <; 2 but a factor 
of two lower (higher) at z ;> 3. These differences are 
consequences of the fact that Galactic-type (SMC-type) 
dust selectively absorbs less (more) ultraviolet radiation 
than the LMC-type dust. Therefore, the high-z solutions 
of Cl a (z) are uncertain (by a factor of two) due to the 
uncertainty in the optical opacity of grains, comparable 
to or less than the uncertainty caused by the small-scale 
dumpiness of the dust. The solutions presented here are 
all based on LMC-type dust, because this is an inter- 
mediate case between Galactic and SMC-type dust and 
because, in combination with a Salpeter IMF, it seems 
to account well for all the available data. 

(5) The infrared opacity of grains, responsible for the 
dust emission, is approximated here by a power law, 
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K a (u) oc v n , where the normalization and index (n = 2) 
are fixed by our adopted LMC-type dust (see Fig. 4). 
The normalization has roughly the same effect as r c y 
on the amplitude of J v , while the index has roughly the 
same effect as e on the shape of J v . Since our fits to the 
cosmic infrared background involve adjustment of t c v 
and e, the uncertainties in the infrared opacity of grains 
are partially absorbed by these parameters. We have 
checked that a different infrared opacity (within 50% 
of the LMC-type dust) would have only minor effects 
on our results (but with different values of T c y and e). 
However, if the normalization were much higher (lower) 
than adopted here, it would require much more (less) 
absorption of stellar energy to reproduce the same in- 
frared background, causing the optical background to 
be too blue (red). A value of n = 1.5 would only give a 
marginally acceptable fit to the far-infrared background. 

4.2. Observational Uncertainties 

We have used as input to our models the best deter- 
minations available in 1997 of the comoving density of 
HI and the rest-frame ultraviolet cmissivity, Ohio and 
E uo . Any errors in Ohio and E vo would propagate into 
the corresponding solutions for Q g and f2 s . Since the 
quoted statistical errors in Jlnio and E vo are all 60% or 
less, these uncertainties do not strongly affect our re- 
sults. In addition to the statistical errors, however, it 
is always possible that the input data are plagued by 
unknown systematic errors and/or biases. After com- 
pleting this work, we became aware of new estimates 
of E vo at high redshifts by Steidel et al. (1999) from 
a large ground-based survey of Lyman-brcak galaxies: 
\o%E vo = 19.42 ± 0.07 at z = 3.04 and logE^ = 
19.33 ± 0.10 at z = 4.13, both for A = 1500 A and 
by integrating a Schechter luminosity function down to 
0.1 £» (Steidel 1999, private communication). The new 
estimate of E vo at z « 4 is higher than the previous one 
by Madau et al. (1998) from the HDF by a factor of two. 
Steidel et al. (1999) suggest that the reason for this dis- 
crepancy is that the HDF may be too small to be a fair 
sample of the universe. Another source of uncertainty 
in E vo is the fact that the faint end of the luminosity 
function has not been well established at high redshifts. 

We have repeated all our calculations with the new 
estimates of E vo by Steidel et al. (1999) at z > 3 and 
the entries in Table 2 at lower redshifts. We then read- 
justed the parameters T c y and e in the models to match 
the cosmic far-infrared background measurements. The 
resulting best, minimum, and maximum solutions for 
fl s are plotted in Figure 12, with (r c y,e) = (0.82,0.8), 
(0.98,1.4), and (0.88,0.4), respectively. These solutions 
are not radically different from those shown in Figure 



7b. The largest difference between the best solutions in 
the two cases is only 50% at z <; 4, even though the 
corresponding input values of E vo at z « 4 differ by a 
factor of two. The reason for this is that our solutions 
for tt s are constrained by many observational input data 
and hence the propagation of errors from E vo to (l s is 
significantly reduced (by a factor of two at z « 4) . Given 
these uncertainties in the input data, the global rate of 
star formation is fairly well determined by our models. 
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FIG. 12. — Comoving rate of star formation in our 
models based on the modified emissivities at z > 3 from 
Steidel et al. (1999). The symbols have the same mean- 
ing as in Fig. 7b. 



4.3. Comparisons with Other Results 

The global rates of star formation derived here are 
compared with independent estimates made by other 
techniques in Figure 13. The solid, short-dash, and 
long-dash curves are respectively the best, minimum, 
and maximum solutions for 0*(z) in our models. The 
dotted curves are the predictions of the inflow (upper 
curve) and outflow (lower curve) models of Pei & Fall 
(1995), which are based on absorption-line observations 
of damped Lya systems. The circles at z — 0, 0.2, and 
0.8, come from Ha observations of galaxies in a local 
sample (Gallego et al. 1995) and in the CFRS sam- 
ple (Tresse & Maddox 1998; Glazebrook et al. 1999). 
For consistency, we have adopted the Glazebrook et 
al. (1999) conversion from the Ha luminosity density 
L-Ra (corrected for dust using H/3-to-Ha line ratios) to 
the star formation rate (appropriate for the Salpeter 
IMF with solar metallicity) , equivalent to (f2*/Gyr _1 ) = 
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1.0 x 10" 43 (i H a/erg s" 1 Mpc~ 3 ). A similar conversion 
factor is also found in our models (assuming 0.45 Ha 
photons per Lyman continuum photon for case B re- 
combination and including the metallicity dependence 
of the spectral evolution). The squares at z « 0.2 — 1.0 
come from ultraviolet, optical, near- infrared, and radio 
observations of ISO 15 /im selected galaxies in the CFRS 
survey (Flores et al. 1999). Clearly, our solutions for the 
global rates of star formation agree well with the avail- 
able Ha and mid- infrared surveys at z <^ 1. They are 
also consistent with an estimate of the star formation 
rate at z « 3 by Hughes et al. (1998) from submillime- 
ter observations of the HDF with SCUBA, although this 
depends on only a few detected sources and assumptions 
about their dust temperatures and redshifts (see addi- 
tional results from Barger et al. 1998). The global his- 
tory of star formation derived here is also broadly con- 
sistent with recent estimates of the redshift distribution 
of SCUBA sources (Lilly et al. 1999). 

In summary, we conclude that our solutions for the 



-5.0 



1 








f 

J 


T 1 

;■ / 

7/ 

4/ 

If : 


11 / 


. \ \ 

\-.\ \ 
\ V^X" 


- T// 

VI 

1 






1 , , , , 1 , , , , 





1 2 3 4 5 



z 

FIG. 13. — Comparisons of the comoving rate of star 
formation in our models with estimates from recent Ha, 
mid-infrared, and submillimeter observations. The solid, 
short-dash, and long-dash curves are respectively the 
best, minimum, and maximum solutions derived here. 
The two dotted curves are the predictions of Pei & Fall 
(1995) from models of cosmic chemical evolution. The 
circles at z = 0, 0.2, and 0.875 are from the Ha sur- 
veys by Gallego et al. (1995), Tresse & Maddox (1998), 
and Glazebrook et al. (1999), respectively. The squares 
are from the ISO 15/um survey by Flores et al. (1999). 
The circle with an arrow at z sa 3 is from the SCUBA 
submillimeter observations by Hughes et al. (1998). 



global rate of star formation are robust at z <J 2, but 
still uncertain by a factor of three at z ;> 3. The ro- 
bustness at z ^ 2 is due to the tight constraint imposed 
by the cosmic infrared background on how much radia- 
tion is emitted and absorbed at shorter wavelengths. At 
z J> 3, there is little time to accumulate much radiation 
in the background. Consequently, the observed back- 
ground provides only a weak constraint on the star for- 
mation rates at high redshifts. Our models indicate that 
the rest-frame emissivity of galaxies evolves strongly not 
only in the near-ultraviolet (e.g., 2800 A) but also in the 
far-infrared (e.g., 100 /im). Either one of these emissivi- 
ties, if studied in isolation, would provide only a partial 
tracer of star formation. The most promising way to 
determine the global rates of star formation at z > 3 
will be to combine future near-infrared and submillime- 
ter surveys; the former traces the direct radiation from 
stars at rest-frame near-ultraviolet wavelengths, while 
the latter traces the reradiation by dust at rest-frame 
far-infrared wavelengths. For z J> 4, stellar emission 
at A — 2800 A is redshifted to A ^ 1.4 fim, and falls 
within the window of the Next Generation Space Tele- 
scope (NGST), whereas dust emission at A = 100 /mi is 
redshifted to A ;> 500 fim, and falls within the window 
of the Millimeter Array (MMA). Thus, a combination 
of these two instruments should eventually enable us to 
trace the global history of star formation to high red- 
shifts. 

4.4. Evolution of Galaxies 

Our models of cosmic chemical evolution provide a 
global description of the major constituents of galaxies. 
Figure 14 shows the inferred evolution of the stars, in- 
terstellar gas, baryons in galaxies, and baryon flow be- 
tween the interstellar and intergalactic media. The data 
point f2 s = (4.9^2^) x 10~ 3 at z = is an estimate 
from the observed mean B-band emissivity of nearby 
galaxies and the mean mass-to-light ratios in spheroids, 
disks, and irregular galaxies (Fukugita, Hogan, & Pee- 
bles 1998), while the data point tt g = (5.0 ± 1.2) x 10~ 4 
at z — comes from observations of the 21 cm emis- 
sion of nearby galaxies (§ 2.2.6). The results from our 
models are plotted as the curves in each panel of Figure 
14. The agreement with (0) is by design, while the 
agreement with n s (0) is an indication that our adopted 
lower mass cutoff in the IMF is reasonable. Evidently, 
galaxies contain mostly stars at z <J 1 and mostly gas 
at z ^ 2. The bulk of the stars seen today formed at 
relatively low redshifts: 34 - 48% since z = 1, 72 - 91% 
since z = 2, 91 - 98% since z = 3, and 97 - 99% since 
z = 4. 

Although there are still significant uncertainties in 
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FIG. 14. — Evolution of (a) the comoving density of stars, (b) the comoving density of interstellar gas, (c) the 
comoving density of baryons in galaxies, and (d) the comoving rate of baryon flow between the interstellar and 
intergalactic media. The solid, short-dash, and long-dash curves are the best, minimum, and maximum solutions in 
our models. The data point for f2 s (0) is an estimate from the mean blue luminosity density and mass-to- light ratio 
of nearby galaxies (Fukugita ct al. 1998), while the data point for f2 s (0) is inferred from 21 cm emission observations 
of nearby galaxies (§ 2.2.6). 



these solutions, they can be divided roughly into three 
periods: (1) At z ;> 3, there is a rapid inflow from the 
intergalactic medium and hence a build up of the in- 
terstellar gas within galaxies (Figs. 14b and 14c). (2) 
At 1 <; z <; 3, the interstellar gas content is high 
(Fig. 14b), presumably triggering high rates of star for- 
mation (Fig. 13). This in turn causes the rapid accumu- 



lation of stars (Fig. 14a) and heavy elements (Fig. 8), 
and strong emission from stars and dust (Fig. 9). (3) 
At z <J 1, the stellar content continues to increase but 
at a slower rate, while the gas content (Fig. 14b), star 
formation rate (Fig. 13), ultraviolet and optical stellar 
emission (Fig. 9), and far- infrared dust emission (Fig. 9) 
all decrease. During this period, the interstellar gas is 
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consumed by star formation without much replenish- 
ment by inflow. Our results therefore suggest that, on 
average, galaxies spend 13 — 18% of their lives (defined 
by the present age of the universe) in a "growth" pe- 
riod, 23 — 25% of their lives in a "working" period, and 
64 — 57% of their lives in a "retirement" period (for A = 
and 0.1 < % < 0.5). 

It appears from our models that the most significant 
inflow occurred at z w 3 — 4 (Fig. 14d). This infcrrence, 
however, depends crucially on the estimates of ftmo from 
damped Lya surveys at the highest redshifts (Fig. 7a). 
If the observed HI content were underestimated at z ;> 
3, the peak in the inflow could be pushed to an earlier 
epoch. Since the interstellar metallicity is inversely pro- 
portional to the comoving density of gas, the observed 
metallicity at z « 3 — 4 also provide a constraint on the 
peak in the inflow. Therefore, better determinations of 
the comoving density of HI and mean metallicity in the 
damped Lya systems at z ;> 3 would help reduce the un- 
certainties in the epoch when the bulk of the interstellar 
gas was assembled in galaxies. 

In most models based on hierarchical clustering, in- 
cluding the cold dark matter model, galaxies form rel- 
atively late. The exact evolution depends on a large 
number of processes, including the merging of dark mat- 
ter halos, the gravitational heating and radiative cooling 
of gas within them, the formation of stars, the feed- 
back of stellar mass, energy, and metals to the interstel- 
lar medium, and so forth. Semianalytical models treat 
these processes by simple, physically-motivated recipes; 
and with suitable choices of parameters, they are found 
to be consistent with many of the observed properties 
of galaxies (White & Frenk 1991; Kauffmann, White, 
Guiderdoni 1993; Cole et al. 1994; Somerville & Pri- 
mack 1998). In some semianalytical models, the comov- 
ing rate of star formation and the comoving density of 
cold gas rise to peak values at z « 1 — 2 and then decline 
at lower redshifts (Kauffmann 1996; Baugh et al. 1998). 
These models are therefore qualitatively consistent with 
our models of cosmic chemical evolution and the obser- 
vations on which they are based. The evolution of ft s 
and fig in the semianalytical models is, however, milder 
than the observed evolution. It is not yet clear whether 
or not this represents a significant discrepancy. 



write the comoving density of metals in the intergalactic 
medium as 
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can be positive or 
We assume, as a good approximation, that 
f for Q f < (by ignoring any inflow if the net 



where fto is the comoving rate of outflow from galaxies 
and Z is the mean metallicity of the outflowing gas. The 
net flow rate discussed in § 2.1 is simply ttf = fii — 
fto, where fti is the comoving rate of inflow. Both Oo 
and fti must be positive, whereas ft/ 
negative 
fto = -ft 

flow is dominated by outflow) and fio = t]fl f for ft/ > 
with ?)Cl (to allow for some outflow if the net flow is 
dominated by inflow) . A small but non-zero value of r\ 
only affects equation (3) in § 2.1. We have checked that 
all of the results presented in § 3 are nearly the same 
for r\ <J 0.3, with a maximum difference of less than 20% 
in Z. In this case, we can make rough estimates of the 
chemical enrichment of the intergalactic medium from 
our solutions for the net rate of baryon flow. 

Figure 15 shows the mean comoving density of metals 
in the intergalactic medium, resulting from the outflow 
of metal-enriched gas from galaxies in our models. The 
three curves represent our three solutions with r\ = 0.2. 
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4.5. Chemical Enrichment of the Intergalactic Medium 

The chemical enrichment of the intergalactic medium 
could be the result of an outflow of metal-enriched gas 
from galaxies or nucleosynthesis in objects much smaller 
and more uniformly distributed than galaxies (e.g., pre- 
galactic stars or star clusters) . We explore in this section 
the possibility of outflows from galaxies. To do so, we 



FIG. 15. — Mean comoving density of metals expelled 
from galaxies into the intergalactic medium as a func- 
tion of rcdshift. The solid, short-dash, and long-dash 
curves are the best, minimum, and maximum solutions 
in our models. The data point at z = 3.2 was derived 
by Songaila (1997) from observations of C IV and Si IC 
absorption lines in the Lya forest with some ionization 
corrections and assumptions about relative abundances. 



26 



The data point at z = 3.2 was derived by Songaila 
(1997) from observations of C IV and Si IV absorption 
lines in the Lya forest (corresponding roughly to Nm <; 
10 17 cm -2 ) with some modest ionization corrections and 
assumptions about the relative abundances of different 
elements. Although the data and models are both un- 
certain, the rough agreement indicates that even a small 
outflow of metal-enriched gas from galaxies could ac- 
count for the observed metal content of the Lya forest. 
A comparison of Figure 15 with Figure 8 suggests that 
the comoving density of intergalactic metals is typically 
an order of magnitude smaller than that of interstellar 
metals and that the two are comparable only at z ^ 0.4. 
This is nearly independent of the relative outflow rate 
for rj <^ 0.3. Galaxies probably retained most of their 
metals up until now, but there could be almost as many 
metals in the intergalactic medium as in the interstellar 
media of galaxies at z = 0. 

The global mean mctallicity of the intergalactic 
medium is (Zi GM ) = fimiGM/^giGM, where ^ ff i G M is the 
total gas content in the intergalactic medium. Adopting 
^giGM ~ 0.02 from Big Bang nucleosynthesis (Fukugita 
ct al. 1998), we estimate log(ZiGM/Z Q ) ~ —1-5, —2.6, 
and —3.6 at z = 0, 2, and 4, respectively. We em- 
phasize, however, that there could be a large dispersion 
around the mean, depending on whether and how much 
the metals expelled from galaxies are mixed with the 
intergalactic medium. If the metals were mixed com- 
pletely with all of the intergalactic gas, the metallicity 
would be the same everywhere, i.e., Ziqm = (-Zigm)- If 
the metals (and gas) were not mixed at all in the in- 
tergalactic medium, the mctallicity would fluctuate be- 
tween Ziqm = Z in regions where the ejected gas is 
present and Ziqm = in regions where it is not present. 
Some degree of inhomogeneity would seem to be the 
most likely situation. 

5. CONCLUSIONS 

We have studied the global evolution of the stellar, 
gaseous, and metal contents of galaxies, as well as the 
radiation they emit and the baryon flow between the in- 
terstellar and intergalactic media. Our results are based 
on models of cosmic chemical evolution, with observa- 
tional inputs from damped Lya absorption-line surveys 
and optical imaging and redshift surveys, and measure- 
ments of the cosmic infrared background. The novel 
feature of our approach is a self-consistent treatment of 
the obscuration, absorption, and emission by dust. We 
emphasize that all of our results pertain to the mean 
properties of galaxies, averaged over comoving volumes 
large enough to be fair samples of galaxies of all masses, 
sizes, and morphological types. Such a global approach 



has the advantage of simplicity, but the disadvantage 
that many interesting questions concerning the diver- 
sity of galaxies cannot be addressed. In particular, our 
results are completely independent of the unknown mor- 
phologies, sizes, and space densities of the damped Lya 
systems. 

We find that the extragalactic far-infrared back- 
ground detected by the COBE DIRBE and FIRAS can 
be accounted for if much of the radiation emitted by 
stars is absorbed and reemitted by small-scale clumps 
of dust, possibly associated with individual stars. The 
ratio of the bolometric dust emission to the bolomet- 
ric stellar emission in our models increases from 36% at 
z = to about 60% at z « 1 and becomes uncertain at 
higher rcdshifts (10 — 70% at z ;> 3). Of the total energy 
radiated by stars (integrated over all redshifts), 41 — 46% 
is in the ultraviolet /optical/near-infrared part of extra- 
galactic background, while the remaining 59 — 54% is in 
the mid-infrarcd/far-infrared/submillimcter part. The 
absorption by dust corrects the global rates of star for- 
mation inferred from optical surveys upward by factors 
of 2.5, 3 — 4, and 2 — 8 at z = 0, 1, and 3, respectively. 
From our models, it appears that galaxies may have ex- 
perienced a sustained period of intense star formation 
at 1 <; z <; 3 and somewhat less rapid star formation at 
z <; 1 and z ;> 3, roughly correlated with the evolution 
of their neutral gas content, although the evolution at 
z ^ 3 is quite uncertain. We note, however, that if the 
global rate of star formation were to remain high beyond 
z « 3, it would produce more metals in the interstellar 
medium than are observed in the damped Lya systems, 
unless a significant outflow had occurred before z ~ 3. 

The solutions presented here for the mean properties 
of galaxies are based on relatively few input observa- 
tions, but they match remarkably well a wide variety of 
other observations, mostly within the quoted lcr errors 
of the data. In particular, our models reproduce: (1) the 
extragalactic background derived from galaxy counts at 
A = 0.2-2.2 ^m; (2) the rest-frame 0.44, 1.0, and 2.2 ^m 
emissivities from the available surveys of galaxies up to 
z w 2; (3) the infrared dust emissivities from nearby 
galaxies at A = 12, 25, 60, and 100 /an; (4) the observed 
mean abundance of heavy elements in the damped Lya 
systems at 0.4 <; z <; 3.5; and (5) the global rates of 
star formation inferred from Ha and ISO mid-infrared 
surveys at z <J 1 and recent SCUBA submillimeter ob- 
servations. The metal enrichment history of the inter- 
galactic medium implied by our models is also consistent 
with the observed mean comoving density of metals in 
the Lyman-alpha forest at high redshifts. Our solutions 
are robust at z <; 2 but uncertain at z ;> 3. Major 
sources of uncertainties are the optical opacity of grains 
and the small-scale dumpiness of dust in the interstellar 
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medium. Both of these are likely to remain uncertain 
for some time and hence will affect the interpretations 
of future optical and infrared observations as well. 

Our results indicate that many different observations 
of galaxies, by means of quasar absorption lines, optical 
imaging and redshift surveys, and extragalactic back- 
ground measurements, can be brought together to pro- 
vide a coherent picture of the global evolution of galaxies 
over much of the Hubble time. The models presented 
here suggest that the process of galaxy formation may 
have undergone different evolutionary phases: (1) an 
early period of significant inflow to assemble interstel- 
lar gas at z ;> 3; (2) a subsequent period of intense star 
formation and chemical enrichment at 1 <J z <J 3; and 
(3) a recent period of decline in the gas content, star 
formation rate, and the ultraviolet and optical stellar 
emissivity, and far- infrared dust emissivity at z <; 1. 
Currently, there are still large uncertainties in this pic- 
ture, especially at z J> 3. Future surveys that aim to 
trace the rest-frame near-ultraviolet emissivity of stars, 
the rest- frame far-infrared emissivity of dust, and the 
global properties of interstellar gas, heavy elements, and 
dust beyond z w 3 could reduce these uncertainties and 
thus advance our knowledge of the global evolution of 
galaxies. 
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